AliPhysics  master (3d17d9d)
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 <TLorentzVector.h>
25 #include "AliMCEvent.h"
26 #include "AliAODEvent.h"
27 #include "AliAODMCHeader.h"
28 #include "AliGenEventHeader.h"
29 #include "AliAODMCParticle.h"
30 #include "AliAODRecoDecayHF.h"
31 #include "AliVertexingHFUtils.h"
32 
33 /* $Id$ */
34 
36 // //
37 // Class with functions useful for different D2H analyses //
38 // - event plane resolution //
39 // - <pt> calculation with side band subtraction //
40 // - tracklet multiplicity calculation //
41 // Origin: F.Prino, Torino, prino@to.infn.it //
42 // //
44 
46 ClassImp(AliVertexingHFUtils);
48 
49 using namespace std;
50 
51 //______________________________________________________________________
53  fK(1),
54  fSubRes(1.),
55  fMinEtaForTracklets(-1.),
56  fMaxEtaForTracklets(1.)
57 {
59 }
60 
61 
62 //______________________________________________________________________
64  TObject(),
65  fK(k),
66  fSubRes(1.),
69 {
71 }
72 
73 
74 //______________________________________________________________________
75 void AliVertexingHFUtils::ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance,Double_t &errsignificance){
77 
78 
79  Double_t errSigSq=errsignal*errsignal;
80  Double_t errBkgSq=errbackground*errbackground;
81  Double_t sigPlusBkg=signal+background;
82  if(sigPlusBkg>0. && signal>0.){
83  significance = signal/TMath::Sqrt(signal+background);
84  errsignificance = significance*TMath::Sqrt((errSigSq+errBkgSq)/(4.*sigPlusBkg*sigPlusBkg)+(background/sigPlusBkg)*errSigSq/signal/signal);
85  }else{
86  significance=0.;
87  errsignificance=0.;
88  }
89  return;
90 
91 }
92 //______________________________________________________________________
95  Double_t c[5];
96  if(k==1){
97  c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
98  }
99  else if(k==2){
100  c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
101  }
102  else return -1;
103  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;
104 }
105 
106 //______________________________________________________________________
108  return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
109 }
110 
111 
112 //______________________________________________________________________
115 
116  Double_t x1=0;
117  Double_t x2=15;
118  Double_t y1,y2;
119  if(k==1){
120  y1=ResolK1(x1)-res;
121  y2=ResolK1(x2)-res;
122  }
123  else if(k==2){
124  y1=Pol(x1,2)-res;
125  y2=Pol(x2,2)-res;
126  }
127  else return -1;
128 
129  if(y1*y2>0) return -1;
130  if(y1==0) return y1;
131  if(y2==0) return y2;
132  Double_t xmed,ymed;
133  Int_t jiter=0;
134  while((x2-x1)>0.0001){
135  xmed=0.5*(x1+x2);
136  if(k==1){
137  y1=ResolK1(x1)-res;
138  y2=ResolK1(x2)-res;
139  ymed=ResolK1(xmed)-res;
140  }
141  else if(k==2){
142  y1=Pol(x1,2)-res;
143  y2=Pol(x2,2)-res;
144  ymed=Pol(xmed,2)-res;
145  }
146  else return -1;
147  if((y1*ymed)<0) x2=xmed;
148  if((y2*ymed)<0)x1=xmed;
149  if(ymed==0) return xmed;
150  jiter++;
151  }
152  return 0.5*(x1+x2);
153 }
154 
155 //______________________________________________________________________
158  Double_t chisub=FindChi(resSub,k);
159  Double_t chifull=chisub*TMath::Sqrt(2);
160  if(k==1) return ResolK1(chifull);
161  else if(k==2) return Pol(chifull,2);
162  else{
163  printf("k should be <=2\n");
164  return 1.;
165  }
166 }
167 
168 //______________________________________________________________________
171  if(!hSubEvCorr) return 1.;
172  Double_t resSub=GetSubEvResol(hSubEvCorr);
173  return GetFullEvResol(resSub,k);
174 }
175 //______________________________________________________________________
178  if(!hSubEvCorr) return 1.;
179  Double_t resSub=GetSubEvResolLowLim(hSubEvCorr);
180  printf("%f\n",resSub);
181  return GetFullEvResol(resSub,k);
182 }
183 //______________________________________________________________________
186  if(!hSubEvCorr) return 1.;
187  Double_t resSub=GetSubEvResolHighLim(hSubEvCorr);
188  printf("%f\n",resSub);
189  return GetFullEvResol(resSub,k);
190 }
191 //______________________________________________________________________
194  AliAODTracklets* tracklets=ev->GetTracklets();
195  Int_t nTr=tracklets->GetNumberOfTracklets();
196  Int_t count=0;
197  for(Int_t iTr=0; iTr<nTr; iTr++){
198  Double_t theta=tracklets->GetTheta(iTr);
199  Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
200  if(eta>mineta && eta<maxeta) count++;
201  }
202  return count;
203 }
204 //______________________________________________________________________
207 
208  Int_t nChargedMC=0;
209  for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
210  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
211  Int_t charge = part->Charge();
212  Double_t eta = part->Eta();
213  if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++;
214  }
215  return nChargedMC;
216 }
217 //______________________________________________________________________
220 
221  Int_t nChargedMC=0;
222  for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
223  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
224  Int_t charge = part->Charge();
225  Double_t eta = part->Eta();
226  if(charge!=0 && eta>mineta && eta<maxeta){
227  if(part->IsPrimary())nChargedMC++;
228  }
229  }
230  return nChargedMC;
231 }
232 //______________________________________________________________________
235 
236  Int_t nChargedMC=0;
237  for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
238  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
239  Int_t charge = part->Charge();
240  Double_t eta = part->Eta();
241  if(charge!=0 && eta>mineta && eta<maxeta){
242  if(part->IsPhysicalPrimary())nChargedMC++;
243  }
244  }
245  return nChargedMC;
246 }
247 
248 
249 //______________________________________________________________________
251  //
254  // http://git.cern.ch/pubweb/AliRoot.git/blob/HEAD:/ANALYSIS/AliCentralitySelectionTask.cxx#l1345
255 
256  Double_t multV0AEq=0;
257  for(Int_t iCh = 32; iCh < 64; ++iCh) {
258  Double_t mult = ev->GetVZEROEqMultiplicity(iCh);
259  multV0AEq += mult;
260  }
261  return multV0AEq;
262 }
263 
264 //______________________________________________________________________
269 
270  Double_t multV0CEq=0;
271  for(Int_t iCh = 0; iCh < 32; ++iCh) {
272  Double_t mult = ev->GetVZEROEqMultiplicity(iCh);
273  multV0CEq += mult;
274  }
275  return multV0CEq;
276 }
277 
278 //______________________________________________________________________
279 void AliVertexingHFUtils::AveragePt(Float_t& averagePt, Float_t& errorPt,Float_t ptmin,Float_t ptmax, TH2F* hMassD, Float_t massFromFit, Float_t sigmaFromFit,
280  TF1* funcB2, Float_t sigmaRangeForSig,Float_t sigmaRangeForBkg, Float_t minMass, Float_t maxMass, Int_t rebin){
281 
283 
284  //Make 2D histos in the desired pt range
285  Int_t start=hMassD->FindBin(ptmin);
286  Int_t end=hMassD->FindBin(ptmax)-1;
287  const Int_t nx=end-start;
288  TH2F *hMassDpt=new TH2F("hptmass","hptmass",nx,ptmin,ptmax,hMassD->GetNbinsY(),hMassD->GetYaxis()->GetBinLowEdge(1),
289  hMassD->GetYaxis()->GetBinLowEdge(hMassD->GetNbinsY())+hMassD->GetYaxis()->GetBinWidth(hMassD->GetNbinsY()));
290  for(Int_t ix=start;ix<end;ix++){
291  for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){
292  hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy));
293  hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy));
294  }
295  }
296 
297  Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
298  Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
299  Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig);
300  Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig);
301  Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig);
302  Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig);
303  // printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
304 
305  Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
306  Int_t minBinBkgLow=TMath::Max(hMassD->GetYaxis()->FindBin(minMass),2);
307  Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
308  Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
309  Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
310  Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
311  Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
312  Int_t maxBinBkgHi=TMath::Min(hMassD->GetYaxis()->FindBin(maxMass),hMassD->GetNbinsY()-1);
313  Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
314  Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
315  // printf("BKG Fit Limits = %f %f && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
316 
317  Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin);
318  Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
319  Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
320  // printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
321 
322  TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow);
323  TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi);
324  TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig);
325 
326  hMptBkgLo->Rebin(rebin);
327  hMptBkgHi->Rebin(rebin);
328  hMptSigReg->Rebin(rebin);
329 
330  hMptBkgLo->Sumw2();
331  hMptBkgHi->Sumw2();
332  TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
333  hMptBkgLoScal->Scale(bkgSig/bkgLow);
334  TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
335  hMptBkgHiScal->Scale(bkgSig/bkgHi);
336 
337  TH1F* hMptBkgAver=0x0;
338  hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin");
339  hMptBkgAver->Add(hMptBkgHiScal);
340  hMptBkgAver->Scale(0.5);
341  TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin");
342  hMptSig->Add(hMptBkgAver,-1.);
343 
344  averagePt = hMptSig->GetMean();
345  errorPt = hMptSig->GetMeanError();
346 
347  delete hMptBkgLo;
348  delete hMptBkgHi;
349  delete hMptSigReg;
350  delete hMptBkgLoScal;
351  delete hMptBkgHiScal;
352  delete hMptBkgAver;
353  delete hMassDpt;
354  delete hMptSig;
355 
356 }
357 //____________________________________________________________________________
360  const Double32_t *mean = aodEv->GetT0TOF();
361  if(mean && mean[0]<9999.) return kTRUE;
362  else return kFALSE;
363 }
364 //____________________________________________________________________________
365 Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
367 
368  if(!partD || !arrayMC || !mcHeader) return 99999.;
369  Int_t code=partD->GetPdgCode();
370  if(TMath::Abs(code)!=421) return 99999.;
371 
372  Double_t vtxTrue[3];
373  mcHeader->GetVertex(vtxTrue);
374  Double_t origD[3];
375  partD->XvYvZv(origD);
376  Short_t charge=partD->Charge();
377  Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
378  for(Int_t iDau=0; iDau<2; iDau++){
379  pXdauTrue[iDau]=0.;
380  pYdauTrue[iDau]=0.;
381  pZdauTrue[iDau]=0.;
382  }
383 
384  Int_t nDau=partD->GetNDaughters();
385  Int_t labelFirstDau = partD->GetDaughterLabel(0);
386  if(nDau==2){
387  for(Int_t iDau=0; iDau<2; iDau++){
388  Int_t ind = labelFirstDau+iDau;
389  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
390  if(!part){
391  printf("Daughter particle not found in MC array");
392  return 99999.;
393  }
394  pXdauTrue[iDau]=part->Px();
395  pYdauTrue[iDau]=part->Py();
396  pZdauTrue[iDau]=part->Pz();
397  }
398  }else{
399  return 99999.;
400  }
401 
402  Double_t d0dummy[2]={0.,0.};
403  AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
404  return aodDvsMC.ImpParXY();
405 
406 }
407 //____________________________________________________________________________
408 Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
410 
411  if(!partD || !arrayMC || !mcHeader) return 99999.;
412  Int_t code=partD->GetPdgCode();
413  if(TMath::Abs(code)!=411) return 99999.;
414 
415  Double_t vtxTrue[3];
416  mcHeader->GetVertex(vtxTrue);
417  Double_t origD[3];
418  partD->XvYvZv(origD);
419  Short_t charge=partD->Charge();
420  Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
421  for(Int_t iDau=0; iDau<3; iDau++){
422  pXdauTrue[iDau]=0.;
423  pYdauTrue[iDau]=0.;
424  pZdauTrue[iDau]=0.;
425  }
426 
427  Int_t nDau=partD->GetNDaughters();
428  Int_t labelFirstDau = partD->GetDaughterLabel(0);
429  if(nDau==3){
430  for(Int_t iDau=0; iDau<3; iDau++){
431  Int_t ind = labelFirstDau+iDau;
432  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
433  if(!part){
434  printf("Daughter particle not found in MC array");
435  return 99999.;
436  }
437  pXdauTrue[iDau]=part->Px();
438  pYdauTrue[iDau]=part->Py();
439  pZdauTrue[iDau]=part->Pz();
440  }
441  }else if(nDau==2){
442  Int_t theDau=0;
443  for(Int_t iDau=0; iDau<2; iDau++){
444  Int_t ind = labelFirstDau+iDau;
445  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
446  if(!part){
447  printf("Daughter particle not found in MC array");
448  return 99999.;
449  }
450  Int_t pdgCode=TMath::Abs(part->GetPdgCode());
451  if(pdgCode==211 || pdgCode==321){
452  pXdauTrue[theDau]=part->Px();
453  pYdauTrue[theDau]=part->Py();
454  pZdauTrue[theDau]=part->Pz();
455  ++theDau;
456  }else{
457  Int_t nDauRes=part->GetNDaughters();
458  if(nDauRes==2){
459  Int_t labelFirstDauRes = part->GetDaughterLabel(0);
460  for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
461  Int_t indDR = labelFirstDauRes+iDauRes;
462  AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
463  if(!partDR){
464  printf("Daughter particle not found in MC array");
465  return 99999.;
466  }
467 
468  Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
469  if(pdgCodeDR==211 || pdgCodeDR==321){
470  pXdauTrue[theDau]=partDR->Px();
471  pYdauTrue[theDau]=partDR->Py();
472  pZdauTrue[theDau]=partDR->Pz();
473  ++theDau;
474  }
475  }
476  }
477  }
478  }
479  if(theDau!=3){
480  printf("Wrong number of decay prongs");
481  return 99999.;
482  }
483  }
484 
485  Double_t d0dummy[3]={0.,0.,0.};
486  AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
487  return aodDvsMC.ImpParXY();
488 
489 }
490 
491 
492 
493 //____________________________________________________________________________
494 Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) {
495  //
496  // Correct the number of accepted tracklets based on a TProfile Hist
497  //
498  //
499 
500  if(TMath::Abs(vtxZ)>10.0){
501  // printf("ERROR: Z vertex out of range for correction of multiplicity\n");
502  return uncorrectedNacc;
503  }
504 
505  if(!estimatorAvg){
506  printf("ERROR: Missing TProfile for correction of multiplicity\n");
507  return uncorrectedNacc;
508  }
509 
510  Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
511 
512  Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
513 
514  Double_t correctedNacc = std::max(0., uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM)));
515 
516  return correctedNacc;
517 }
518 //______________________________________________________________________
519 TString AliVertexingHFUtils::GetGenerator(Int_t label, AliAODMCHeader* header){
521 
522  Int_t nsumpart=0;
523  TList *lh=header->GetCocktailHeaders();
524  Int_t nh=lh->GetEntries();
525  for(Int_t i=0;i<nh;i++){
526  AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
527  TString genname=gh->GetName();
528  Int_t npart=gh->NProduced();
529  if(label>=nsumpart && label<(nsumpart+npart)) return genname;
530  nsumpart+=npart;
531  }
532  TString empty="";
533  return empty;
534 }
535 //_____________________________________________________________________
536 void AliVertexingHFUtils::GetTrackPrimaryGenerator(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
537  GetTrackPrimaryGenerator(track->GetLabel(),header,arrayMC,nameGen);
538 }
539 //_____________________________________________________________________
540 void AliVertexingHFUtils::GetTrackPrimaryGenerator(Int_t label, AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
541 
543 
544  Int_t lab=TMath::Abs(label);
545  nameGen=GetGenerator(lab,header);
546 
547  // Int_t countControl=0;
548 
549  while(nameGen.IsWhitespace()){
550  AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
551  if(!mcpart) break;
552  Int_t mother = mcpart->GetMother();
553  if(mother<0) break;
554  lab=mother;
555  nameGen=GetGenerator(mother,header);
556  // countControl++;
557  // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
558  // printf("AliVertexingHFUtils::GetTrackPrimaryGenerator - BREAK: Protection from infinite loop active\n");
559  // break;
560  // }
561  }
562 
563  return;
564 }
565 //----------------------------------------------------------------------
566 Bool_t AliVertexingHFUtils::IsTrackInjected(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC){
568 
569  return IsTrackInjected(track->GetLabel(),header,arrayMC);
570 }
571 //----------------------------------------------------------------------
572 Bool_t AliVertexingHFUtils::IsTrackInjected(Int_t label, AliAODMCHeader *header,TClonesArray *arrayMC){
574  TString nameGen;
575  Int_t lab=TMath::Abs(label);
576 
577  GetTrackPrimaryGenerator(lab,header,arrayMC,nameGen);
578 
579  if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
580 
581  return kTRUE;
582 }
583 //____________________________________________________________________________
584 Bool_t AliVertexingHFUtils::IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
587  Int_t nprongs=cand->GetNProngs();
588  for(Int_t i=0;i<nprongs;i++){
589  AliAODTrack *daugh=(AliAODTrack*)cand->GetDaughter(i);
590  if(IsTrackInjected(daugh,header,arrayMC)) return kTRUE;
591  }
592  return kFALSE;
593 }
594 //____________________________________________________________________________
595 Bool_t AliVertexingHFUtils::IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODEvent* aod, AliAODMCHeader *header,TClonesArray *arrayMC){
599 
600  Int_t nprongs=cand->GetNProngs();
601  for(Int_t i=0;i<nprongs;i++){
602  Int_t idDau=cand->GetProngID(i);
603  for(Int_t i=0; i<aod->GetNumberOfTracks(); i++) {
604  AliAODTrack *daugh=(AliAODTrack*)aod->GetTrack(i);
605  if(daugh && (Int_t)daugh->GetID()==idDau){
606  if(AliVertexingHFUtils::IsTrackInjected(daugh,header,arrayMC)) return kTRUE;
607  }
608  }
609  }
610  return kFALSE;
611 }
612 //____________________________________________________________________________
613 Bool_t AliVertexingHFUtils::HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
615 
616  AliAODTrack* bach = cand->GetBachelor();
617  if(IsTrackInjected(bach, header, arrayMC)) {
618  // printf("Bachelor is injected, the whole candidate is then injected\n");
619  return kTRUE;
620  }
621  AliAODv0* v0 = cand->Getv0();
622  Int_t nprongs = v0->GetNProngs();
623  for(Int_t i = 0; i < nprongs; i++){
624  AliAODTrack *daugh = (AliAODTrack*)v0->GetDaughter(i);
625  if(IsTrackInjected(daugh,header,arrayMC)) {
626  // printf("V0 daughter number %d is injected, the whole candidate is then injected\n", i);
627  return kTRUE;
628  }
629  }
630  return kFALSE;
631 }
632 //____________________________________________________________________________
633 Int_t AliVertexingHFUtils::PreSelectITSUpgrade(TClonesArray* arrayMC, AliAODMCHeader *header, TObjArray aodTracks, Int_t nDaug, Int_t pdgabs, const Int_t *pdgDg){
639 
640  Bool_t injected = kFALSE;
641  Bool_t hijing = kFALSE;
642 
643  AliAODTrack *track[(const Int_t)nDaug];
644  Int_t dgLabels[(const Int_t)nDaug];
645  for(Int_t iD=0; iD<nDaug; iD++) {
646  track[iD] = (AliAODTrack*)aodTracks.At(iD);
647  dgLabels[iD] = track[iD]->GetLabel();
648 
649  Bool_t isTrInjected = IsTrackInjected(track[iD],header,arrayMC);
650  if(isTrInjected) injected = kTRUE;
651  else hijing = kTRUE;
652 
653  //Combination of injected signal + HIJING background
654  if(injected && hijing) return 0;
655  }
656 
657  //Purely HIJING background, no need to MatchToMC
658  if(hijing) return 3;
659 
660  //Purely injected signal, check by Matching to MC if it is considered as signal
661  // Mimic MatchToMC (can't access as it is protected), code below is the same as:
662  // Int_t MatchedToLabel = AliAODRecoDecay::MatchToMC(pdgabs, arrayMC, dgLabels, nDaug, 0, pdgDg);
663  if(injected){
664  //Label of mother particle, or -1 when some operation failed.
665  Int_t MatchedToLabel = 0;
666 
667  Int_t labMom[10]={0,0,0,0,0,0,0,0,0,0};
668  Int_t i,j,lab,labMother,pdgMother,pdgPart;
669  AliAODMCParticle *part=0;
670  AliAODMCParticle *mother=0;
671  Bool_t pdgUsed[10]={kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
672 
673  // loop on daughter labels
674  for(i=0; i<nDaug; i++) {
675  if(MatchedToLabel == -1) continue;
676 
677  labMom[i]=-1;
678  lab = TMath::Abs(dgLabels[i]);
679  if(lab<0){ MatchedToLabel = -1; continue; }
680 
681  part = (AliAODMCParticle*)arrayMC->At(lab);
682  if(!part){ MatchedToLabel = -1; continue; }
683 
684  // check the PDG of the daughter, if requested
685  pdgPart=TMath::Abs(part->GetPdgCode());
686  for(j=0; j<nDaug; j++) {
687  if(!pdgUsed[j] && pdgPart==pdgDg[j]) {
688  pdgUsed[j]=kTRUE;
689  break;
690  }
691  }
692 
693  mother = part;
694  while(mother->GetMother()>=0) {
695  labMother=mother->GetMother();
696  mother = (AliAODMCParticle*)arrayMC->At(labMother);
697  if(!mother) break;
698 
699  pdgMother = TMath::Abs(mother->GetPdgCode());
700  if(pdgMother==pdgabs) {
701  labMom[i]=labMother;
702  break;
703  } else if(pdgMother>pdgabs || pdgMother<10) {
704  break;
705  }
706  }
707  if(labMom[i]==-1) MatchedToLabel = -1; // mother PDG not ok for this daughter
708  } // end loop on daughters
709 
710  if(MatchedToLabel == 0){
711  // check if the candidate is signal
712  labMother=labMom[0];
713  // all labels have to be the same and !=-1
714  for(i=0; i<nDaug; i++) {
715  if(labMom[i]==-1) MatchedToLabel = -1;
716  if(labMom[i]!=labMother) MatchedToLabel = -1;
717  }
718 
719  // check that all daughter PDGs are matched
720  for(i=0; i<nDaug; i++) {
721  if(pdgUsed[i]==kFALSE) MatchedToLabel = -1;
722  }
723  }
724 
725  if(MatchedToLabel==0) MatchedToLabel = labMother;
726 
727  if(MatchedToLabel != -1) return 1; //injected, matched to MC
728  else return 2; //injected, not matched to MC
729  }
730 
731  //Should not reach this, if so check. Return 0 for compilation warning
732  printf("AliVertexingHFUtils::PreSelectITSUpgrade: Neither injected, nor HIJING");
733  return 0;
734 }
735 //____________________________________________________________________________
736 Int_t AliVertexingHFUtils::CheckOrigin(AliMCEvent* mcEvent, AliMCParticle *mcPart, Bool_t searchUpToQuark){
738 
739  Int_t pdgGranma = 0;
740  Int_t mother = 0;
741  mother = mcPart->GetMother();
742  Int_t istep = 0;
743  Int_t abspdgGranma =0;
744  Bool_t isFromB=kFALSE;
745  Bool_t isQuarkFound=kFALSE;
746  while (mother >=0 ){
747  istep++;
748  AliMCParticle* mcGranma = (AliMCParticle*)mcEvent->GetTrack(mother);
749  if (mcGranma){
750  TParticle* partGranma = mcGranma->Particle();
751  if(partGranma) pdgGranma = partGranma->GetPdgCode();
752  abspdgGranma = TMath::Abs(pdgGranma);
753  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
754  isFromB=kTRUE;
755  }
756  if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
757  mother = mcGranma->GetMother();
758  }else{
759  printf("CheckOrigin: Failed casting the mother particle!");
760  break;
761  }
762  }
763  if(searchUpToQuark && !isQuarkFound) return 0;
764  if(isFromB) return 5;
765  else return 4;
766 
767 }
768 //____________________________________________________________________________
769 Int_t AliVertexingHFUtils::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark){
771 
772  Int_t pdgGranma = 0;
773  Int_t mother = 0;
774  mother = mcPart->GetMother();
775  Int_t istep = 0;
776  Int_t abspdgGranma =0;
777  Bool_t isFromB=kFALSE;
778  Bool_t isQuarkFound=kFALSE;
779  while (mother >=0 ){
780  istep++;
781  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
782  if (mcGranma){
783  pdgGranma = mcGranma->GetPdgCode();
784  abspdgGranma = TMath::Abs(pdgGranma);
785  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
786  isFromB=kTRUE;
787  }
788  if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
789  mother = mcGranma->GetMother();
790  }else{
791  printf("AliVertexingHFUtils::CheckOrigin: Failed casting the mother particle!");
792  break;
793  }
794  }
795  if(searchUpToQuark && !isQuarkFound) return 0;
796  if(isFromB) return 5;
797  else return 4;
798 
799 }
800 //____________________________________________________________________________
801 Bool_t AliVertexingHFUtils::IsTrackFromCharm(AliAODTrack* tr, TClonesArray* arrayMC){
803  Int_t absLabel=TMath::Abs(tr->GetLabel());
804  AliAODMCParticle* mcPart=dynamic_cast<AliAODMCParticle*>(arrayMC->At(absLabel));
805  Int_t mother = mcPart->GetMother();
806  Int_t istep = 0;
807  while (mother >=0 ){
808  istep++;
809  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
810  if (mcGranma){
811  Int_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
812  if ((abspdgGranma==4) ||(abspdgGranma>400 && abspdgGranma<500) || (abspdgGranma>4000 && abspdgGranma<5000)) return kTRUE;
813  mother = mcGranma->GetMother();
814  }else{
815  printf("AliVertexingHFUtils::IsTrackFromCharm: Failed casting the mother particle!");
816  break;
817  }
818  }
819  return kFALSE;
820 }
821 //____________________________________________________________________________
822 Bool_t AliVertexingHFUtils::IsTrackFromBeauty(AliAODTrack* tr, TClonesArray* arrayMC){
824  Int_t absLabel=TMath::Abs(tr->GetLabel());
825  AliAODMCParticle* mcPart=dynamic_cast<AliAODMCParticle*>(arrayMC->At(absLabel));
826  Int_t mother = mcPart->GetMother();
827  Int_t istep = 0;
828  while (mother >=0 ){
829  istep++;
830  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
831  if (mcGranma){
832  Int_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
833  if ((abspdgGranma==5) ||(abspdgGranma>500 && abspdgGranma<600) || (abspdgGranma>5000 && abspdgGranma<6000)) return kTRUE;
834  mother = mcGranma->GetMother();
835  }else{
836  printf("AliVertexingHFUtils::IsTrackFromBeauty: Failed casting the mother particle!");
837  break;
838  }
839  }
840  return kFALSE;
841 }
842 //____________________________________________________________________________
843 Bool_t AliVertexingHFUtils::IsTrackFromHadronDecay(Int_t pdgMoth, AliAODTrack* tr, TClonesArray* arrayMC){
845  Int_t absLabel=TMath::Abs(tr->GetLabel());
846  Int_t absPdgMoth=TMath::Abs(pdgMoth);
847  AliAODMCParticle* mcPart=dynamic_cast<AliAODMCParticle*>(arrayMC->At(absLabel));
848  Int_t mother = mcPart->GetMother();
849  Int_t istep = 0;
850  while (mother >=0 ){
851  istep++;
852  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
853  if (mcGranma){
854  Int_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
855  if (abspdgGranma==absPdgMoth) return kTRUE;
856  mother = mcGranma->GetMother();
857  }else{
858  printf("AliVertexingHFUtils::IsTrackFromHadronDecay: Failed casting the mother particle!");
859  break;
860  }
861  }
862  return kFALSE;
863 }
864 //____________________________________________________________________________
865 Double_t AliVertexingHFUtils::GetBeautyMotherPt(TClonesArray* arrayMC, AliAODMCParticle *mcPart){
867 
868  Int_t pdgGranma = 0;
869  Int_t mother = 0;
870  mother = mcPart->GetMother();
871  Int_t istep = 0;
872  Int_t abspdgGranma =0;
873  while (mother >=0 ){
874  istep++;
875  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
876  if (mcGranma){
877  pdgGranma = mcGranma->GetPdgCode();
878  abspdgGranma = TMath::Abs(pdgGranma);
879  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
880  return mcGranma->Pt();
881  }
882  if(abspdgGranma==4) return -999.;
883  if(abspdgGranma==5) return -1.;
884  mother = mcGranma->GetMother();
885  }else{
886  printf("AliVertexingHFUtils::GetBeautyMotherPt: Failed casting the mother particle!");
887  break;
888  }
889  }
890  return -999.;
891 }
892 
893 //____________________________________________________________________________
894 Double_t AliVertexingHFUtils::GetBeautyMotherPtAndPDG(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t &pdgGranma){
896 
897  pdgGranma = 0;
898  Int_t mother = 0;
899  mother = mcPart->GetMother();
900  Int_t istep = 0;
901  Int_t abspdgGranma = 0;
902  while (mother >=0 ){
903  istep++;
904  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
905  if (mcGranma){
906  pdgGranma = mcGranma->GetPdgCode();
907  abspdgGranma = TMath::Abs(pdgGranma);
908  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
909  return mcGranma->Pt();
910  }
911  if(abspdgGranma==4) return -999.;
912  if(abspdgGranma==5) return -1.;
913  mother = mcGranma->GetMother();
914  }else{
915  printf("AliVertexingHFUtils::GetBeautyMotherPt: Failed casting the mother particle!");
916  break;
917  }
918  }
919  pdgGranma = 0;
920  return -999.;
921 }
922 
923 //____________________________________________________________________________
924 Int_t AliVertexingHFUtils::CheckD0Decay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
926 
927  if(label<0) return -1;
928  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
929  TParticle* part = mcEvent->Particle(label);
930  if(!part || !mcPart) return -1;
931  Int_t pdgD=part->GetPdgCode();
932  if(TMath::Abs(pdgD)!=421) return -1;
933 
934  Int_t nDau=mcPart->GetNDaughters();
935 
936  if(nDau==2){
937  Int_t daughter0 = mcPart->GetDaughterFirst();
938  Int_t daughter1 = mcPart->GetDaughterLast();
939  TParticle* mcPartDaughter0 = mcEvent->Particle(daughter0);
940  TParticle* mcPartDaughter1 = mcEvent->Particle(daughter1);
941  if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
942  arrayDauLab[0]=daughter0;
943  arrayDauLab[1]=daughter1;
944  Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
945  Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
946  if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
947  !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
948  return -1;
949  }
950  if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
951  if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
952  if((pdgCode0*pdgCode1)>0) return -1;
953  Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
954  Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
955  Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
956  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
957  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
958  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
959  return 1;
960  }
961 
962  if(nDau==3 || nDau==4){
963  Int_t nKaons=0;
964  Int_t nPions=0;
965  Double_t sumPxDau=0.;
966  Double_t sumPyDau=0.;
967  Double_t sumPzDau=0.;
968  Int_t nFoundKpi=0;
969  Int_t labelFirstDau = mcPart->GetDaughterFirst();
970  for(Int_t iDau=0; iDau<nDau; iDau++){
971  Int_t indDau = labelFirstDau+iDau;
972  if(indDau<0) return -1;
973  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
974  TParticle* dau=mcEvent->Particle(indDau);
975  if(!mcDau || !dau) return -1;
976  Int_t pdgdau=dau->GetPdgCode();
977  if(TMath::Abs(pdgdau)==321){
978  if(pdgD>0 && pdgdau>0) return -1;
979  if(pdgD<0 && pdgdau<0) return -1;
980  nKaons++;
981  sumPxDau+=dau->Px();
982  sumPyDau+=dau->Py();
983  sumPzDau+=dau->Pz();
984  arrayDauLab[nFoundKpi++]=indDau;
985  if(nFoundKpi>4) return -1;
986  }else if(TMath::Abs(pdgdau)==211){
987  nPions++;
988  sumPxDau+=dau->Px();
989  sumPyDau+=dau->Py();
990  sumPzDau+=dau->Pz();
991  arrayDauLab[nFoundKpi++]=indDau;
992  if(nFoundKpi>4) return -1;
993  }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
994  Int_t nResDau=mcDau->GetNDaughters();
995  if(nResDau!=2) return -1;
996  Int_t indFirstResDau=mcDau->GetDaughterFirst();
997  for(Int_t resDau=0; resDau<2; resDau++){
998  Int_t indResDau=indFirstResDau+resDau;
999  if(indResDau<0) return -1;
1000  TParticle* resdau=mcEvent->Particle(indResDau);
1001  if(!resdau) return -1;
1002  Int_t pdgresdau=resdau->GetPdgCode();
1003  if(TMath::Abs(pdgresdau)==321){
1004  if(pdgD>0 && pdgresdau>0) return -1;
1005  if(pdgD<0 && pdgresdau<0) return -1;
1006  nKaons++;
1007  sumPxDau+=resdau->Px();
1008  sumPyDau+=resdau->Py();
1009  sumPzDau+=resdau->Pz();
1010  arrayDauLab[nFoundKpi++]=indResDau;
1011  if(nFoundKpi>4) return -1;
1012  }
1013  if(TMath::Abs(pdgresdau)==211){
1014  nPions++;
1015  sumPxDau+=resdau->Px();
1016  sumPyDau+=resdau->Py();
1017  sumPzDau+=resdau->Pz();
1018  arrayDauLab[nFoundKpi++]=indResDau;
1019  if(nFoundKpi>4) return -1;
1020  }
1021  }
1022  }else{
1023  return -1;
1024  }
1025  }
1026  if(nPions!=3) return -1;
1027  if(nKaons!=1) return -1;
1028  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -1;
1029  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -1;
1030  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -1;
1031  return 2;
1032  }
1033  return -1;
1034 }
1035 
1036 //____________________________________________________________________________
1037 Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1039 
1040  Int_t pdgD=mcPart->GetPdgCode();
1041  if(TMath::Abs(pdgD)!=421) return -1;
1042 
1043  Int_t nDau=mcPart->GetNDaughters();
1044 
1045  if(nDau==2){
1046  Int_t daughter0 = mcPart->GetDaughterLabel(0);
1047  Int_t daughter1 = mcPart->GetDaughterLabel(1);
1048  AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter0));
1049  AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter1));
1050  if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
1051  arrayDauLab[0]=daughter0;
1052  arrayDauLab[1]=daughter1;
1053  Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
1054  Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
1055  if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
1056  !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
1057  return -1;
1058  }
1059  if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
1060  if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
1061  if((pdgCode0*pdgCode1)>0) return -1;
1062  Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
1063  Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
1064  Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
1065  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1066  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1067  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1068  return 1;
1069  }
1070 
1071  if(nDau==3 || nDau==4){
1072  Int_t nKaons=0;
1073  Int_t nPions=0;
1074  Double_t sumPxDau=0.;
1075  Double_t sumPyDau=0.;
1076  Double_t sumPzDau=0.;
1077  Int_t nFoundKpi=0;
1078  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
1079  for(Int_t iDau=0; iDau<nDau; iDau++){
1080  Int_t indDau = labelFirstDau+iDau;
1081  if(indDau<0) return -1;
1082  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1083  if(!dau) return -1;
1084  Int_t pdgdau=dau->GetPdgCode();
1085  if(TMath::Abs(pdgdau)==321){
1086  if(pdgD>0 && pdgdau>0) return -1;
1087  if(pdgD<0 && pdgdau<0) return -1;
1088  nKaons++;
1089  sumPxDau+=dau->Px();
1090  sumPyDau+=dau->Py();
1091  sumPzDau+=dau->Pz();
1092  arrayDauLab[nFoundKpi++]=indDau;
1093  if(nFoundKpi>4) return -1;
1094  }else if(TMath::Abs(pdgdau)==211){
1095  nPions++;
1096  sumPxDau+=dau->Px();
1097  sumPyDau+=dau->Py();
1098  sumPzDau+=dau->Pz();
1099  arrayDauLab[nFoundKpi++]=indDau;
1100  if(nFoundKpi>4) return -1;
1101  }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
1102  Int_t nResDau=dau->GetNDaughters();
1103  if(nResDau!=2) return -1;
1104  Int_t indFirstResDau=dau->GetDaughterLabel(0);
1105  for(Int_t resDau=0; resDau<2; resDau++){
1106  Int_t indResDau=indFirstResDau+resDau;
1107  if(indResDau<0) return -1;
1108  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1109  if(!resdau) return -1;
1110  Int_t pdgresdau=resdau->GetPdgCode();
1111  if(TMath::Abs(pdgresdau)==321){
1112  if(pdgD>0 && pdgresdau>0) return -1;
1113  if(pdgD<0 && pdgresdau<0) return -1;
1114  nKaons++;
1115  sumPxDau+=resdau->Px();
1116  sumPyDau+=resdau->Py();
1117  sumPzDau+=resdau->Pz();
1118  arrayDauLab[nFoundKpi++]=indResDau;
1119  if(nFoundKpi>4) return -1;
1120  }
1121  if(TMath::Abs(pdgresdau)==211){
1122  nPions++;
1123  sumPxDau+=resdau->Px();
1124  sumPyDau+=resdau->Py();
1125  sumPzDau+=resdau->Pz();
1126  arrayDauLab[nFoundKpi++]=indResDau;
1127  if(nFoundKpi>4) return -1;
1128  }
1129  }
1130  }else{
1131  return -1;
1132  }
1133  }
1134  if(nPions!=3) return -1;
1135  if(nKaons!=1) return -1;
1136  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
1137  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
1138  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
1139  return 2;
1140  }
1141  return -1;
1142 }
1143 //____________________________________________________________________________
1144 Int_t AliVertexingHFUtils::CheckDplusDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1146 
1147  if(label<0) return -1;
1148  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1149  TParticle* part = mcEvent->Particle(label);
1150  if(!part || !mcPart) return -1;
1151  Int_t pdgD=part->GetPdgCode();
1152  if(TMath::Abs(pdgD)!=411) return -1;
1153 
1154  Int_t nDau=mcPart->GetNDaughters();
1155  Int_t labelFirstDau = mcPart->GetDaughterFirst();
1156  Int_t nKaons=0;
1157  Int_t nPions=0;
1158  Double_t sumPxDau=0.;
1159  Double_t sumPyDau=0.;
1160  Double_t sumPzDau=0.;
1161  Int_t nFoundKpi=0;
1162 
1163  if(nDau==3 || nDau==2){
1164  for(Int_t iDau=0; iDau<nDau; iDau++){
1165  Int_t indDau = labelFirstDau+iDau;
1166  if(indDau<0) return -1;
1167  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1168  TParticle* dau=mcEvent->Particle(indDau);
1169  if(!mcDau || !dau) return -1;
1170  Int_t pdgdau=dau->GetPdgCode();
1171  if(TMath::Abs(pdgdau)==321){
1172  if(pdgD*pdgdau>0) return -1;
1173  nKaons++;
1174  sumPxDau+=dau->Px();
1175  sumPyDau+=dau->Py();
1176  sumPzDau+=dau->Pz();
1177  arrayDauLab[nFoundKpi++]=indDau;
1178  if(nFoundKpi>3) return -1;
1179  }else if(TMath::Abs(pdgdau)==211){
1180  if(pdgD*pdgdau<0) return -1;
1181  nPions++;
1182  sumPxDau+=dau->Px();
1183  sumPyDau+=dau->Py();
1184  sumPzDau+=dau->Pz();
1185  arrayDauLab[nFoundKpi++]=indDau;
1186  if(nFoundKpi>3) return -1;
1187  }else if(TMath::Abs(pdgdau)==313){
1188  Int_t nResDau=mcDau->GetNDaughters();
1189  if(nResDau!=2) return -1;
1190  Int_t indFirstResDau=mcDau->GetDaughterFirst();
1191  for(Int_t resDau=0; resDau<2; resDau++){
1192  Int_t indResDau=indFirstResDau+resDau;
1193  if(indResDau<0) return -1;
1194  TParticle* resdau=mcEvent->Particle(indResDau);
1195  if(!resdau) return -1;
1196  Int_t pdgresdau=resdau->GetPdgCode();
1197  if(TMath::Abs(pdgresdau)==321){
1198  if(pdgD*pdgresdau>0) return -1;
1199  sumPxDau+=resdau->Px();
1200  sumPyDau+=resdau->Py();
1201  sumPzDau+=resdau->Pz();
1202  nKaons++;
1203  arrayDauLab[nFoundKpi++]=indResDau;
1204  if(nFoundKpi>3) return -1;
1205  }
1206  if(TMath::Abs(pdgresdau)==211){
1207  if(pdgD*pdgresdau<0) return -1;
1208  sumPxDau+=resdau->Px();
1209  sumPyDau+=resdau->Py();
1210  sumPzDau+=resdau->Pz();
1211  nPions++;
1212  arrayDauLab[nFoundKpi++]=indResDau;
1213  if(nFoundKpi>3) return -1;
1214  }
1215  }
1216  }else{
1217  return -1;
1218  }
1219  }
1220  if(nPions!=2) return -1;
1221  if(nKaons!=1) return -1;
1222  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1223  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1224  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1225  if(nDau==3) return 1;
1226  else if(nDau==2) return 2;
1227  }
1228 
1229  return -1;
1230 
1231 }
1232 
1233 //____________________________________________________________________________
1234 Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1236 
1237  Int_t pdgD=mcPart->GetPdgCode();
1238  if(TMath::Abs(pdgD)!=411) return -1;
1239 
1240  Int_t nDau=mcPart->GetNDaughters();
1241  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
1242  Int_t nKaons=0;
1243  Int_t nPions=0;
1244  Double_t sumPxDau=0.;
1245  Double_t sumPyDau=0.;
1246  Double_t sumPzDau=0.;
1247  Int_t nFoundKpi=0;
1248 
1249  if(nDau==3 || nDau==2){
1250  for(Int_t iDau=0; iDau<nDau; iDau++){
1251  Int_t indDau = labelFirstDau+iDau;
1252  if(indDau<0) return -1;
1253  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1254  if(!dau) return -1;
1255  Int_t pdgdau=dau->GetPdgCode();
1256  if(TMath::Abs(pdgdau)==321){
1257  if(pdgD*pdgdau>0) return -1;
1258  nKaons++;
1259  sumPxDau+=dau->Px();
1260  sumPyDau+=dau->Py();
1261  sumPzDau+=dau->Pz();
1262  arrayDauLab[nFoundKpi++]=indDau;
1263  if(nFoundKpi>3) return -1;
1264  }else if(TMath::Abs(pdgdau)==211){
1265  if(pdgD*pdgdau<0) return -1;
1266  nPions++;
1267  sumPxDau+=dau->Px();
1268  sumPyDau+=dau->Py();
1269  sumPzDau+=dau->Pz();
1270  arrayDauLab[nFoundKpi++]=indDau;
1271  if(nFoundKpi>3) return -1;
1272  }else if(TMath::Abs(pdgdau)==313){
1273  Int_t nResDau=dau->GetNDaughters();
1274  if(nResDau!=2) return -1;
1275  Int_t indFirstResDau=dau->GetDaughterLabel(0);
1276  for(Int_t resDau=0; resDau<2; resDau++){
1277  Int_t indResDau=indFirstResDau+resDau;
1278  if(indResDau<0) return -1;
1279  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1280  if(!resdau) return -1;
1281  Int_t pdgresdau=resdau->GetPdgCode();
1282  if(TMath::Abs(pdgresdau)==321){
1283  if(pdgD*pdgresdau>0) return -1;
1284  sumPxDau+=resdau->Px();
1285  sumPyDau+=resdau->Py();
1286  sumPzDau+=resdau->Pz();
1287  nKaons++;
1288  arrayDauLab[nFoundKpi++]=indResDau;
1289  if(nFoundKpi>3) return -1;
1290  }
1291  if(TMath::Abs(pdgresdau)==211){
1292  if(pdgD*pdgresdau<0) return -1;
1293  sumPxDau+=resdau->Px();
1294  sumPyDau+=resdau->Py();
1295  sumPzDau+=resdau->Pz();
1296  nPions++;
1297  arrayDauLab[nFoundKpi++]=indResDau;
1298  if(nFoundKpi>3) return -1;
1299  }
1300  }
1301  }else{
1302  return -1;
1303  }
1304  }
1305  if(nPions!=2) return -1;
1306  if(nKaons!=1) return -1;
1307  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1308  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1309  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1310  if(nDau==3) return 1;
1311  else if(nDau==2) return 2;
1312  }
1313 
1314  return -1;
1315 
1316 }
1317 //____________________________________________________________________________
1318 Int_t AliVertexingHFUtils::CheckDplusKKpiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1320 
1321  if(label<0) return -1;
1322  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1323  TParticle* part = mcEvent->Particle(label);
1324  if(!part || !mcPart) return -1;
1325  Int_t pdgD=part->GetPdgCode();
1326  if(TMath::Abs(pdgD)!=411) return -1;
1327 
1328  Int_t nDau=mcPart->GetNDaughters();
1329  Int_t labelFirstDau = mcPart->GetDaughterFirst();
1330  Int_t nKaons=0;
1331  Int_t nPions=0;
1332  Double_t sumPxDau=0.;
1333  Double_t sumPyDau=0.;
1334  Double_t sumPzDau=0.;
1335  Int_t nFoundKpi=0;
1336  Bool_t isPhi=kFALSE;
1337  Bool_t isk0st=kFALSE;
1338 
1339  if(nDau==3 || nDau==2){
1340  for(Int_t iDau=0; iDau<nDau; iDau++){
1341  Int_t indDau = labelFirstDau+iDau;
1342  if(indDau<0) return -1;
1343  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1344  TParticle* dau=mcEvent->Particle(indDau);
1345  if(!mcDau || !dau) return -1;
1346  Int_t pdgdau=dau->GetPdgCode();
1347  if(TMath::Abs(pdgdau)==321){
1348  nKaons++;
1349  sumPxDau+=dau->Px();
1350  sumPyDau+=dau->Py();
1351  sumPzDau+=dau->Pz();
1352  arrayDauLab[nFoundKpi++]=indDau;
1353  if(nFoundKpi>3) return -1;
1354  }else if(TMath::Abs(pdgdau)==211){
1355  nPions++;
1356  sumPxDau+=dau->Px();
1357  sumPyDau+=dau->Py();
1358  sumPzDau+=dau->Pz();
1359  arrayDauLab[nFoundKpi++]=indDau;
1360  if(nFoundKpi>3) return -1;
1361  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1362  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1363  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1364  Int_t nResDau=mcDau->GetNDaughters();
1365  if(nResDau!=2) return -1;
1366  Int_t indFirstResDau=mcDau->GetDaughterFirst();
1367  for(Int_t resDau=0; resDau<2; resDau++){
1368  Int_t indResDau=indFirstResDau+resDau;
1369  if(indResDau<0) return -1;
1370  TParticle* resdau=mcEvent->Particle(indResDau);
1371  if(!resdau) return -1;
1372  Int_t pdgresdau=resdau->GetPdgCode();
1373  if(TMath::Abs(pdgresdau)==321){
1374  sumPxDau+=resdau->Px();
1375  sumPyDau+=resdau->Py();
1376  sumPzDau+=resdau->Pz();
1377  nKaons++;
1378  arrayDauLab[nFoundKpi++]=indResDau;
1379  if(nFoundKpi>3) return -1;
1380  }
1381  if(TMath::Abs(pdgresdau)==211){
1382  sumPxDau+=resdau->Px();
1383  sumPyDau+=resdau->Py();
1384  sumPzDau+=resdau->Pz();
1385  nPions++;
1386  arrayDauLab[nFoundKpi++]=indResDau;
1387  if(nFoundKpi>3) return -1;
1388  }
1389  }
1390  }else{
1391  return -1;
1392  }
1393  }
1394  if(nPions!=1) return -1;
1395  if(nKaons!=2) return -1;
1396  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1397  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1398  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1399  if(nDau==3) return 3;
1400  else if(nDau==2){
1401  if(isk0st) return 2;
1402  if(isPhi) return 1;
1403  }
1404  }
1405 
1406  return -1;
1407 
1408 }
1409 //____________________________________________________________________________
1410 Int_t AliVertexingHFUtils::CheckDplusKKpiDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1412 
1413  Int_t pdgD=mcPart->GetPdgCode();
1414  if(TMath::Abs(pdgD)!=411) return -1;
1415 
1416  Int_t nDau=mcPart->GetNDaughters();
1417  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
1418  Int_t nKaons=0;
1419  Int_t nPions=0;
1420  Double_t sumPxDau=0.;
1421  Double_t sumPyDau=0.;
1422  Double_t sumPzDau=0.;
1423  Int_t nFoundKpi=0;
1424  Bool_t isPhi=kFALSE;
1425  Bool_t isk0st=kFALSE;
1426 
1427  if(nDau==3 || nDau==2){
1428  for(Int_t iDau=0; iDau<nDau; iDau++){
1429  Int_t indDau = labelFirstDau+iDau;
1430  if(indDau<0) return -1;
1431  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1432  if(!dau) return -1;
1433  Int_t pdgdau=dau->GetPdgCode();
1434  if(TMath::Abs(pdgdau)==321){
1435  nKaons++;
1436  sumPxDau+=dau->Px();
1437  sumPyDau+=dau->Py();
1438  sumPzDau+=dau->Pz();
1439  arrayDauLab[nFoundKpi++]=indDau;
1440  if(nFoundKpi>3) return -1;
1441  }else if(TMath::Abs(pdgdau)==211){
1442  nPions++;
1443  sumPxDau+=dau->Px();
1444  sumPyDau+=dau->Py();
1445  sumPzDau+=dau->Pz();
1446  arrayDauLab[nFoundKpi++]=indDau;
1447  if(nFoundKpi>3) return -1;
1448  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1449  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1450  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1451  Int_t nResDau=dau->GetNDaughters();
1452  if(nResDau!=2) return -1;
1453  Int_t indFirstResDau=dau->GetDaughterLabel(0);
1454  for(Int_t resDau=0; resDau<2; resDau++){
1455  Int_t indResDau=indFirstResDau+resDau;
1456  if(indResDau<0) return -1;
1457  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1458  if(!resdau) return -1;
1459  Int_t pdgresdau=resdau->GetPdgCode();
1460  if(TMath::Abs(pdgresdau)==321){
1461  sumPxDau+=resdau->Px();
1462  sumPyDau+=resdau->Py();
1463  sumPzDau+=resdau->Pz();
1464  nKaons++;
1465  arrayDauLab[nFoundKpi++]=indResDau;
1466  if(nFoundKpi>3) return -1;
1467  }
1468  if(TMath::Abs(pdgresdau)==211){
1469  sumPxDau+=resdau->Px();
1470  sumPyDau+=resdau->Py();
1471  sumPzDau+=resdau->Pz();
1472  nPions++;
1473  arrayDauLab[nFoundKpi++]=indResDau;
1474  if(nFoundKpi>3) return -1;
1475  }
1476  }
1477  }else{
1478  return -1;
1479  }
1480  }
1481  if(nPions!=1) return -1;
1482  if(nKaons!=2) return -1;
1483  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1484  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1485  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1486  if(nDau==3) return 3;
1487  else if(nDau==2){
1488  if(isk0st) return 2;
1489  if(isPhi) return 1;
1490  }
1491  }
1492 
1493  return -1;
1494 
1495 }
1496 //____________________________________________________________________________
1497 Int_t AliVertexingHFUtils::CheckDplusK0spiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1499 
1500  if(label<0) return -1;
1501  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1502  TParticle* part = mcEvent->Particle(label);
1503  if(!part || !mcPart) return -1;
1504  Int_t pdgD=part->GetPdgCode();
1505  if(TMath::Abs(pdgD)!=411) return -1;
1506 
1507  Int_t nDau=mcPart->GetNDaughters();
1508  Int_t labelFirstDau = mcPart->GetDaughterFirst();
1509  Int_t nPions=0;
1510  Double_t sumPxDau=0.;
1511  Double_t sumPyDau=0.;
1512  Double_t sumPzDau=0.;
1513  Int_t nFoundpi=0;
1514 
1515  Int_t codeV0=-1;
1516  if(nDau==2){
1517  for(Int_t iDau=0; iDau<nDau; iDau++){
1518  Int_t indDau = labelFirstDau+iDau;
1519  if(indDau<0) return -1;
1520  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1521  TParticle* dau=mcEvent->Particle(indDau);
1522  if(!mcDau || !dau) return -1;
1523  Int_t pdgdau=dau->GetPdgCode();
1524  if(TMath::Abs(pdgdau)==211){
1525  nPions++;
1526  sumPxDau+=dau->Px();
1527  sumPyDau+=dau->Py();
1528  sumPzDau+=dau->Pz();
1529  arrayDauLab[nFoundpi++]=indDau;
1530  if(nFoundpi>3) return -1;
1531  }else if(TMath::Abs(pdgdau)==311){
1532  codeV0=TMath::Abs(pdgdau);
1533  TParticle* v0=dau;
1534  AliMCParticle* mcV0=mcDau;
1535  if(codeV0==311){
1536  Int_t nK0Dau=mcDau->GetNDaughters();
1537  if(nK0Dau!=1) return -1;
1538  Int_t indK0s=mcDau->GetDaughterFirst();
1539  if(indK0s<0) return -1;
1540  v0=mcEvent->Particle(indK0s);
1541  mcV0=(AliMCParticle*)mcEvent->GetTrack(indK0s);
1542  if(!v0 || !mcV0) return -1;
1543  Int_t pdgK0sl=v0->GetPdgCode();
1544  codeV0=TMath::Abs(pdgK0sl);
1545  }
1546  Int_t nV0Dau=mcV0->GetNDaughters();
1547  if(nV0Dau!=2) return -1;
1548  Int_t indFirstV0Dau=mcV0->GetDaughterFirst();
1549  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
1550  Int_t indV0Dau=indFirstV0Dau+v0Dau;
1551  if(indV0Dau<0) return -1;
1552  TParticle* v0dau=mcEvent->Particle(indV0Dau);
1553  if(!v0dau) return -1;
1554  Int_t pdgv0dau=v0dau->GetPdgCode();
1555  if(TMath::Abs(pdgv0dau)==211){
1556  sumPxDau+=v0dau->Px();
1557  sumPyDau+=v0dau->Py();
1558  sumPzDau+=v0dau->Pz();
1559  nPions++;
1560  arrayDauLab[nFoundpi++]=indV0Dau;
1561  if(nFoundpi>3) return -1;
1562  }
1563  }
1564  }else{
1565  return -1;
1566  }
1567  }
1568  if(nPions!=3) return -1;
1569  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1570  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1571  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1572  if(codeV0==310) return 1;
1573  }
1574  return -1;
1575 
1576 }
1577 
1578 
1579 //____________________________________________________________________________
1580 Int_t AliVertexingHFUtils::CheckDsDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1582 
1583  if(label<0) return -1;
1584  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1585  TParticle* part = mcEvent->Particle(label);
1586  if(!part || !mcPart) return -1;
1587  Int_t pdgD=part->GetPdgCode();
1588  if(TMath::Abs(pdgD)!=431) return -1;
1589 
1590  Int_t nDau=mcPart->GetNDaughters();
1591  Int_t labelFirstDau = mcPart->GetDaughterFirst();
1592  Int_t nKaons=0;
1593  Int_t nPions=0;
1594  Double_t sumPxDau=0.;
1595  Double_t sumPyDau=0.;
1596  Double_t sumPzDau=0.;
1597  Int_t nFoundKpi=0;
1598  Bool_t isPhi=kFALSE;
1599  Bool_t isk0st=kFALSE;
1600  Bool_t isf0=kFALSE;
1601 
1602  if(nDau==3 || nDau==2){
1603  for(Int_t iDau=0; iDau<nDau; iDau++){
1604  Int_t indDau = labelFirstDau+iDau;
1605  if(indDau<0) return -1;
1606  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1607  TParticle* dau=mcEvent->Particle(indDau);
1608  if(!mcDau || !dau) return -1;
1609  Int_t pdgdau=dau->GetPdgCode();
1610  if(TMath::Abs(pdgdau)==321){
1611  nKaons++;
1612  sumPxDau+=dau->Px();
1613  sumPyDau+=dau->Py();
1614  sumPzDau+=dau->Pz();
1615  arrayDauLab[nFoundKpi++]=indDau;
1616  if(nFoundKpi>3) return -1;
1617  }else if(TMath::Abs(pdgdau)==211){
1618  nPions++;
1619  sumPxDau+=dau->Px();
1620  sumPyDau+=dau->Py();
1621  sumPzDau+=dau->Pz();
1622  arrayDauLab[nFoundKpi++]=indDau;
1623  if(nFoundKpi>3) return -1;
1624  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333 || TMath::Abs(pdgdau)==9010221){
1625  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1626  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1627  if(TMath::Abs(pdgdau)==9010221) isf0=kTRUE;
1628  Int_t nResDau=mcDau->GetNDaughters();
1629  if(nResDau!=2) return -1;
1630  Int_t indFirstResDau=mcDau->GetDaughterFirst();
1631  for(Int_t resDau=0; resDau<2; resDau++){
1632  Int_t indResDau=indFirstResDau+resDau;
1633  if(indResDau<0) return -1;
1634  TParticle* resdau=mcEvent->Particle(indResDau);
1635  if(!resdau) return -1;
1636  Int_t pdgresdau=resdau->GetPdgCode();
1637  if(TMath::Abs(pdgresdau)==321){
1638  sumPxDau+=resdau->Px();
1639  sumPyDau+=resdau->Py();
1640  sumPzDau+=resdau->Pz();
1641  nKaons++;
1642  arrayDauLab[nFoundKpi++]=indResDau;
1643  if(nFoundKpi>3) return -1;
1644  }
1645  if(TMath::Abs(pdgresdau)==211){
1646  sumPxDau+=resdau->Px();
1647  sumPyDau+=resdau->Py();
1648  sumPzDau+=resdau->Pz();
1649  nPions++;
1650  arrayDauLab[nFoundKpi++]=indResDau;
1651  if(nFoundKpi>3) return -1;
1652  }
1653  }
1654  }else{
1655  return -1;
1656  }
1657  }
1658  if(nPions!=1) return -1;
1659  if(nKaons!=2) return -1;
1660  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1661  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1662  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1663  if(nDau==3) return 3;
1664  else if(nDau==2){
1665  if(isk0st) return 2;
1666  if(isPhi) return 1;
1667  if(isf0) return 4;
1668  }
1669  }
1670 
1671  return -1;
1672 }
1673 //____________________________________________________________________________
1674 Int_t AliVertexingHFUtils::CheckDsDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1676 
1677  Int_t pdgD=mcPart->GetPdgCode();
1678  if(TMath::Abs(pdgD)!=431) return -1;
1679 
1680  Int_t nDau=mcPart->GetNDaughters();
1681  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
1682  Int_t nKaons=0;
1683  Int_t nPions=0;
1684  Double_t sumPxDau=0.;
1685  Double_t sumPyDau=0.;
1686  Double_t sumPzDau=0.;
1687  Int_t nFoundKpi=0;
1688  Bool_t isPhi=kFALSE;
1689  Bool_t isk0st=kFALSE;
1690  Bool_t isf0=kFALSE;
1691 
1692  if(nDau==3 || nDau==2){
1693  for(Int_t iDau=0; iDau<nDau; iDau++){
1694  Int_t indDau = labelFirstDau+iDau;
1695  if(indDau<0) return -1;
1696  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1697  if(!dau) return -1;
1698  Int_t pdgdau=dau->GetPdgCode();
1699  if(TMath::Abs(pdgdau)==321){
1700  nKaons++;
1701  sumPxDau+=dau->Px();
1702  sumPyDau+=dau->Py();
1703  sumPzDau+=dau->Pz();
1704  arrayDauLab[nFoundKpi++]=indDau;
1705  if(nFoundKpi>3) return -1;
1706  }else if(TMath::Abs(pdgdau)==211){
1707  nPions++;
1708  sumPxDau+=dau->Px();
1709  sumPyDau+=dau->Py();
1710  sumPzDau+=dau->Pz();
1711  arrayDauLab[nFoundKpi++]=indDau;
1712  if(nFoundKpi>3) return -1;
1713  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333 || TMath::Abs(pdgdau)==9010221){
1714  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1715  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1716  if(TMath::Abs(pdgdau)==9010221) isf0=kTRUE;
1717  Int_t nResDau=dau->GetNDaughters();
1718  if(nResDau!=2) return -1;
1719  Int_t indFirstResDau=dau->GetDaughterLabel(0);
1720  for(Int_t resDau=0; resDau<2; resDau++){
1721  Int_t indResDau=indFirstResDau+resDau;
1722  if(indResDau<0) return -1;
1723  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1724  if(!resdau) return -1;
1725  Int_t pdgresdau=resdau->GetPdgCode();
1726  if(TMath::Abs(pdgresdau)==321){
1727  sumPxDau+=resdau->Px();
1728  sumPyDau+=resdau->Py();
1729  sumPzDau+=resdau->Pz();
1730  nKaons++;
1731  arrayDauLab[nFoundKpi++]=indResDau;
1732  if(nFoundKpi>3) return -1;
1733  }
1734  if(TMath::Abs(pdgresdau)==211){
1735  sumPxDau+=resdau->Px();
1736  sumPyDau+=resdau->Py();
1737  sumPzDau+=resdau->Pz();
1738  nPions++;
1739  arrayDauLab[nFoundKpi++]=indResDau;
1740  if(nFoundKpi>3) return -1;
1741  }
1742  }
1743  }else{
1744  return -1;
1745  }
1746  }
1747  if(nPions!=1) return -1;
1748  if(nKaons!=2) return -1;
1749  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1750  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1751  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1752  if(nDau==3) return 3;
1753  else if(nDau==2){
1754  if(isk0st) return 2;
1755  if(isPhi) return 1;
1756  if(isf0) return 4;
1757  }
1758  }
1759 
1760  return -1;
1761 }
1762 //____________________________________________________________________________
1763 Int_t AliVertexingHFUtils::CheckDsK0sKDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1765 
1766  if(label<0) return -1;
1767  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1768  TParticle* part = mcEvent->Particle(label);
1769  if(!part || !mcPart) return -1;
1770  Int_t pdgD=part->GetPdgCode();
1771  if(TMath::Abs(pdgD)!=431) return -1;
1772 
1773  Int_t nDau=mcPart->GetNDaughters();
1774  Int_t labelFirstDau = mcPart->GetDaughterFirst();
1775  Int_t nPions=0;
1776  Int_t nKaons=0;
1777  Double_t sumPxDau=0.;
1778  Double_t sumPyDau=0.;
1779  Double_t sumPzDau=0.;
1780  Int_t nFoundKpi=0;
1781 
1782  Int_t codeV0=-1;
1783  if(nDau==2){
1784  for(Int_t iDau=0; iDau<nDau; iDau++){
1785  Int_t indDau = labelFirstDau+iDau;
1786  if(indDau<0) return -1;
1787  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1788  TParticle* dau=mcEvent->Particle(indDau);
1789  if(!mcDau || !dau) return -1;
1790  Int_t pdgdau=dau->GetPdgCode();
1791  if(TMath::Abs(pdgdau)==211){
1792  nPions++;
1793  sumPxDau+=dau->Px();
1794  sumPyDau+=dau->Py();
1795  sumPzDau+=dau->Pz();
1796  arrayDauLab[nFoundKpi++]=indDau;
1797  if(nFoundKpi>3) return -1;
1798  }else if(TMath::Abs(pdgdau)==321){
1799  nKaons++;
1800  sumPxDau+=dau->Px();
1801  sumPyDau+=dau->Py();
1802  sumPzDau+=dau->Pz();
1803  arrayDauLab[nFoundKpi++]=indDau;
1804  if(nFoundKpi>3) return -1;
1805  }else if(TMath::Abs(pdgdau)==311){
1806  codeV0=TMath::Abs(pdgdau);
1807  TParticle* v0=dau;
1808  AliMCParticle* mcV0=mcDau;
1809  if(codeV0==311){
1810  Int_t nK0Dau=mcDau->GetNDaughters();
1811  if(nK0Dau!=1) return -1;
1812  Int_t indK0s=mcDau->GetDaughterFirst();
1813  if(indK0s<0) return -1;
1814  v0=mcEvent->Particle(indK0s);
1815  mcV0=(AliMCParticle*)mcEvent->GetTrack(indK0s);
1816  if(!v0 || !mcV0) return -1;
1817  Int_t pdgK0sl=v0->GetPdgCode();
1818  codeV0=TMath::Abs(pdgK0sl);
1819  }
1820  Int_t nV0Dau=mcV0->GetNDaughters();
1821  if(nV0Dau!=2) return -1;
1822  Int_t indFirstV0Dau=mcV0->GetDaughterFirst();
1823  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
1824  Int_t indV0Dau=indFirstV0Dau+v0Dau;
1825  if(indV0Dau<0) return -1;
1826  TParticle* v0dau=mcEvent->Particle(indV0Dau);
1827  if(!v0dau) return -1;
1828  Int_t pdgv0dau=v0dau->GetPdgCode();
1829  if(TMath::Abs(pdgv0dau)==211){
1830  sumPxDau+=v0dau->Px();
1831  sumPyDau+=v0dau->Py();
1832  sumPzDau+=v0dau->Pz();
1833  nPions++;
1834  arrayDauLab[nFoundKpi++]=indV0Dau;
1835  if(nFoundKpi>3) return -1;
1836  }
1837  }
1838  }else{
1839  return -1;
1840  }
1841  }
1842  if(nPions!=2) return -1;
1843  if(nKaons!=1) return -1;
1844  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1845  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1846  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1847  if(codeV0==310) return 1;
1848  }
1849  return -1;
1850 
1851 }
1852 //____________________________________________________________________________
1853 Int_t AliVertexingHFUtils::CheckDstarDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1855 
1856  if(label<0) return -1;
1857  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1858  TParticle* part = mcEvent->Particle(label);
1859  if(!part || !mcPart) return -1;
1860  Int_t pdgD=part->GetPdgCode();
1861  if(TMath::Abs(pdgD)!=413) return -1;
1862 
1863  Int_t nDau=mcPart->GetNDaughters();
1864  if(nDau!=2) return -1;
1865 
1866  Int_t labelFirstDau = mcPart->GetDaughterFirst();
1867  Int_t nKaons=0;
1868  Int_t nPions=0;
1869  Double_t sumPxDau=0.;
1870  Double_t sumPyDau=0.;
1871  Double_t sumPzDau=0.;
1872  Int_t nFoundKpi=0;
1873 
1874  for(Int_t iDau=0; iDau<nDau; iDau++){
1875  Int_t indDau = labelFirstDau+iDau;
1876  if(indDau<0) return -1;
1877  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1878  TParticle* dau=mcEvent->Particle(indDau);
1879  if(!mcDau || !dau) return -1;
1880  Int_t pdgdau=dau->GetPdgCode();
1881  if(TMath::Abs(pdgdau)==421){
1882  Int_t nResDau=mcDau->GetNDaughters();
1883  if(nResDau!=2) return -1;
1884  Int_t indFirstResDau=mcDau->GetDaughterFirst();
1885  for(Int_t resDau=0; resDau<2; resDau++){
1886  Int_t indResDau=indFirstResDau+resDau;
1887  if(indResDau<0) return -1;
1888  TParticle* resdau=mcEvent->Particle(indResDau);
1889  if(!resdau) return -1;
1890  Int_t pdgresdau=resdau->GetPdgCode();
1891  if(TMath::Abs(pdgresdau)==321){
1892  if(pdgD*pdgresdau>0) return -1;
1893  sumPxDau+=resdau->Px();
1894  sumPyDau+=resdau->Py();
1895  sumPzDau+=resdau->Pz();
1896  nKaons++;
1897  arrayDauLab[nFoundKpi++]=indResDau;
1898  if(nFoundKpi>3) return -1;
1899  }
1900  if(TMath::Abs(pdgresdau)==211){
1901  if(pdgD*pdgresdau<0) return -1;
1902  sumPxDau+=resdau->Px();
1903  sumPyDau+=resdau->Py();
1904  sumPzDau+=resdau->Pz();
1905  nPions++;
1906  arrayDauLab[nFoundKpi++]=indResDau;
1907  if(nFoundKpi>3) return -1;
1908  }
1909  }
1910  }else if(TMath::Abs(pdgdau)==211){
1911  if(pdgD*pdgdau<0) return -1;
1912  nPions++;
1913  sumPxDau+=dau->Px();
1914  sumPyDau+=dau->Py();
1915  sumPzDau+=dau->Pz();
1916  arrayDauLab[nFoundKpi++]=indDau;
1917  if(nFoundKpi>3) return -1;
1918  }
1919  }
1920 
1921  if(nPions!=2) return -1;
1922  if(nKaons!=1) return -1;
1923  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1924  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1925  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1926  return 1;
1927 
1928 }
1929 
1930 //____________________________________________________________________________
1931 Int_t AliVertexingHFUtils::CheckDstarDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1933 
1934  Int_t pdgD=mcPart->GetPdgCode();
1935  if(TMath::Abs(pdgD)!=413) return -1;
1936 
1937  Int_t nDau=mcPart->GetNDaughters();
1938  if(nDau!=2) return -1;
1939 
1940  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
1941  Int_t nKaons=0;
1942  Int_t nPions=0;
1943  Double_t sumPxDau=0.;
1944  Double_t sumPyDau=0.;
1945  Double_t sumPzDau=0.;
1946  Int_t nFoundKpi=0;
1947 
1948  for(Int_t iDau=0; iDau<nDau; iDau++){
1949  Int_t indDau = labelFirstDau+iDau;
1950  if(indDau<0) return -1;
1951  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1952  if(!dau) return -1;
1953  Int_t pdgdau=dau->GetPdgCode();
1954  if(TMath::Abs(pdgdau)==421){
1955  Int_t nResDau=dau->GetNDaughters();
1956  if(nResDau!=2) return -1;
1957  Int_t indFirstResDau=dau->GetDaughterLabel(0);
1958  for(Int_t resDau=0; resDau<2; resDau++){
1959  Int_t indResDau=indFirstResDau+resDau;
1960  if(indResDau<0) return -1;
1961  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1962  if(!resdau) return -1;
1963  Int_t pdgresdau=resdau->GetPdgCode();
1964  if(TMath::Abs(pdgresdau)==321){
1965  if(pdgD*pdgresdau>0) return -1;
1966  sumPxDau+=resdau->Px();
1967  sumPyDau+=resdau->Py();
1968  sumPzDau+=resdau->Pz();
1969  nKaons++;
1970  arrayDauLab[nFoundKpi++]=indResDau;
1971  if(nFoundKpi>3) return -1;
1972  }
1973  if(TMath::Abs(pdgresdau)==211){
1974  if(pdgD*pdgresdau<0) return -1;
1975  sumPxDau+=resdau->Px();
1976  sumPyDau+=resdau->Py();
1977  sumPzDau+=resdau->Pz();
1978  nPions++;
1979  arrayDauLab[nFoundKpi++]=indResDau;
1980  if(nFoundKpi>3) return -1;
1981  }
1982  }
1983  }else if(TMath::Abs(pdgdau)==211){
1984  if(pdgD*pdgdau<0) return -1;
1985  nPions++;
1986  sumPxDau+=dau->Px();
1987  sumPyDau+=dau->Py();
1988  sumPzDau+=dau->Pz();
1989  arrayDauLab[nFoundKpi++]=indDau;
1990  if(nFoundKpi>3) return -1;
1991  }
1992  }
1993 
1994  if(nPions!=2) return -1;
1995  if(nKaons!=1) return -1;
1996  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1997  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1998  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1999  return 1;
2000 
2001 }
2002 
2003 //____________________________________________________________________________
2004 Int_t AliVertexingHFUtils::CheckLcpKpiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
2006 
2007  if(label<0) return -1;
2008  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
2009  TParticle* part = mcEvent->Particle(label);
2010  if(!part || !mcPart) return -1;
2011  Int_t pdgD=part->GetPdgCode();
2012  if(TMath::Abs(pdgD)!=4122) return -1;
2013 
2014  Int_t nDau=mcPart->GetNDaughters();
2015  Int_t labelFirstDau = mcPart->GetDaughterFirst();
2016  Int_t nKaons=0;
2017  Int_t nPions=0;
2018  Int_t nProtons=0;
2019  Double_t sumPxDau=0.;
2020  Double_t sumPyDau=0.;
2021  Double_t sumPzDau=0.;
2022  Int_t nFoundpKpi=0;
2023 
2024  Int_t codeRes=-1;
2025  if(nDau==3 || nDau==2){
2026  for(Int_t iDau=0; iDau<nDau; iDau++){
2027  Int_t indDau = labelFirstDau+iDau;
2028  if(indDau<0) return -1;
2029  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
2030  TParticle* dau=mcEvent->Particle(indDau);
2031  if(!mcDau || !dau) return -1;
2032  Int_t pdgdau=dau->GetPdgCode();
2033  if(TMath::Abs(pdgdau)==321){
2034  nKaons++;
2035  sumPxDau+=dau->Px();
2036  sumPyDau+=dau->Py();
2037  sumPzDau+=dau->Pz();
2038  arrayDauLab[nFoundpKpi++]=indDau;
2039  if(nFoundpKpi>3) return -1;
2040  }else if(TMath::Abs(pdgdau)==211){
2041  nPions++;
2042  sumPxDau+=dau->Px();
2043  sumPyDau+=dau->Py();
2044  sumPzDau+=dau->Pz();
2045  arrayDauLab[nFoundpKpi++]=indDau;
2046  if(nFoundpKpi>3) return -1;
2047  }else if(TMath::Abs(pdgdau)==2212){
2048  nProtons++;
2049  sumPxDau+=dau->Px();
2050  sumPyDau+=dau->Py();
2051  sumPzDau+=dau->Pz();
2052  arrayDauLab[nFoundpKpi++]=indDau;
2053  if(nFoundpKpi>3) return -1;
2054  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 ||
2055  TMath::Abs(pdgdau)==2224){
2056  codeRes=TMath::Abs(pdgdau);
2057  Int_t nResDau=mcDau->GetNDaughters();
2058  if(nResDau!=2) return -1;
2059  Int_t indFirstResDau=mcDau->GetDaughterFirst();
2060  for(Int_t resDau=0; resDau<2; resDau++){
2061  Int_t indResDau=indFirstResDau+resDau;
2062  if(indResDau<0) return -1;
2063  TParticle* resdau=mcEvent->Particle(indResDau);
2064  if(!resdau) return -1;
2065  Int_t pdgresdau=resdau->GetPdgCode();
2066  if(TMath::Abs(pdgresdau)==321){
2067  sumPxDau+=resdau->Px();
2068  sumPyDau+=resdau->Py();
2069  sumPzDau+=resdau->Pz();
2070  nKaons++;
2071  arrayDauLab[nFoundpKpi++]=indResDau;
2072  if(nFoundpKpi>3) return -1;
2073  }else if(TMath::Abs(pdgresdau)==211){
2074  sumPxDau+=resdau->Px();
2075  sumPyDau+=resdau->Py();
2076  sumPzDau+=resdau->Pz();
2077  nPions++;
2078  arrayDauLab[nFoundpKpi++]=indResDau;
2079  if(nFoundpKpi>3) return -1;
2080  }else if(TMath::Abs(pdgresdau)==2212){
2081  sumPxDau+=resdau->Px();
2082  sumPyDau+=resdau->Py();
2083  sumPzDau+=resdau->Pz();
2084  nProtons++;
2085  arrayDauLab[nFoundpKpi++]=indResDau;
2086  if(nFoundpKpi>3) return -1;
2087  }
2088  }
2089  }else{
2090  return -1;
2091  }
2092  }
2093  if(nPions!=1) return -1;
2094  if(nKaons!=1) return -1;
2095  if(nProtons!=1) return -1;
2096  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
2097  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
2098  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
2099  if(nDau==3) return 1;
2100  else if(nDau==2){
2101  if(codeRes==313) return 2;
2102  else if(codeRes==2224) return 3;
2103  else if(codeRes==3124) return 4;
2104  }
2105  }
2106  return -1;
2107 
2108 }
2109 
2110 //____________________________________________________________________________
2111 Int_t AliVertexingHFUtils::CheckLcpKpiDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
2113 
2114  Int_t pdgD=mcPart->GetPdgCode();
2115  if(TMath::Abs(pdgD)!=4122) return -1;
2116 
2117  Int_t nDau=mcPart->GetNDaughters();
2118  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
2119  Int_t nKaons=0;
2120  Int_t nPions=0;
2121  Int_t nProtons=0;
2122  Double_t sumPxDau=0.;
2123  Double_t sumPyDau=0.;
2124  Double_t sumPzDau=0.;
2125  Int_t nFoundpKpi=0;
2126 
2127  Int_t codeRes=-1;
2128  if(nDau==3 || nDau==2){
2129  for(Int_t iDau=0; iDau<nDau; iDau++){
2130  Int_t indDau = labelFirstDau+iDau;
2131  if(indDau<0) return -1;
2132  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
2133  if(!dau) return -1;
2134  Int_t pdgdau=dau->GetPdgCode();
2135  if(TMath::Abs(pdgdau)==321){
2136  nKaons++;
2137  sumPxDau+=dau->Px();
2138  sumPyDau+=dau->Py();
2139  sumPzDau+=dau->Pz();
2140  arrayDauLab[nFoundpKpi++]=indDau;
2141  if(nFoundpKpi>3) return -1;
2142  }else if(TMath::Abs(pdgdau)==211){
2143  nPions++;
2144  sumPxDau+=dau->Px();
2145  sumPyDau+=dau->Py();
2146  sumPzDau+=dau->Pz();
2147  arrayDauLab[nFoundpKpi++]=indDau;
2148  if(nFoundpKpi>3) return -1;
2149  }else if(TMath::Abs(pdgdau)==2212){
2150  nProtons++;
2151  sumPxDau+=dau->Px();
2152  sumPyDau+=dau->Py();
2153  sumPzDau+=dau->Pz();
2154  arrayDauLab[nFoundpKpi++]=indDau;
2155  if(nFoundpKpi>3) return -1;
2156  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 ||
2157  TMath::Abs(pdgdau)==2224){
2158  codeRes=TMath::Abs(pdgdau);
2159  Int_t nResDau=dau->GetNDaughters();
2160  if(nResDau!=2) return -1;
2161  Int_t indFirstResDau=dau->GetDaughterLabel(0);
2162  for(Int_t resDau=0; resDau<2; resDau++){
2163  Int_t indResDau=indFirstResDau+resDau;
2164  if(indResDau<0) return -1;
2165  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
2166  if(!resdau) return -1;
2167  Int_t pdgresdau=resdau->GetPdgCode();
2168  if(TMath::Abs(pdgresdau)==321){
2169  sumPxDau+=resdau->Px();
2170  sumPyDau+=resdau->Py();
2171  sumPzDau+=resdau->Pz();
2172  nKaons++;
2173  arrayDauLab[nFoundpKpi++]=indResDau;
2174  if(nFoundpKpi>3) return -1;
2175  }else if(TMath::Abs(pdgresdau)==211){
2176  sumPxDau+=resdau->Px();
2177  sumPyDau+=resdau->Py();
2178  sumPzDau+=resdau->Pz();
2179  nPions++;
2180  arrayDauLab[nFoundpKpi++]=indResDau;
2181  if(nFoundpKpi>3) return -1;
2182  }else if(TMath::Abs(pdgresdau)==2212){
2183  sumPxDau+=resdau->Px();
2184  sumPyDau+=resdau->Py();
2185  sumPzDau+=resdau->Pz();
2186  nProtons++;
2187  arrayDauLab[nFoundpKpi++]=indResDau;
2188  if(nFoundpKpi>3) return -1;
2189  }
2190  }
2191  }else{
2192  return -1;
2193  }
2194  }
2195  if(nPions!=1) return -1;
2196  if(nKaons!=1) return -1;
2197  if(nProtons!=1) return -1;
2198  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
2199  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
2200  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
2201  if(nDau==3) return 1;
2202  else if(nDau==2){
2203  if(codeRes==313) return 2;
2204  else if(codeRes==2224) return 3;
2205  else if(codeRes==3124) return 4;
2206  }
2207  }
2208  return -1;
2209 
2210 }
2211 //____________________________________________________________________________
2212 Int_t AliVertexingHFUtils::CheckLcV0bachelorDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
2214 
2215  if(label<0) return -1;
2216  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
2217  TParticle* part = mcEvent->Particle(label);
2218  if(!part || !mcPart) return -1;
2219  Int_t pdgD=part->GetPdgCode();
2220  if(TMath::Abs(pdgD)!=4122) return -1;
2221 
2222  Int_t nDau=mcPart->GetNDaughters();
2223  Int_t labelFirstDau = mcPart->GetDaughterFirst();
2224  Int_t nPions=0;
2225  Int_t nProtons=0;
2226  Double_t sumPxDau=0.;
2227  Double_t sumPyDau=0.;
2228  Double_t sumPzDau=0.;
2229  Int_t nFoundppi=0;
2230 
2231  Int_t codeV0=-1;
2232  if(nDau==2){
2233  for(Int_t iDau=0; iDau<nDau; iDau++){
2234  Int_t indDau = labelFirstDau+iDau;
2235  if(indDau<0) return -1;
2236  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
2237  TParticle* dau=mcEvent->Particle(indDau);
2238  if(!mcDau || !dau) return -1;
2239  Int_t pdgdau=dau->GetPdgCode();
2240  if(TMath::Abs(pdgdau)==211){
2241  nPions++;
2242  sumPxDau+=dau->Px();
2243  sumPyDau+=dau->Py();
2244  sumPzDau+=dau->Pz();
2245  arrayDauLab[nFoundppi++]=indDau;
2246  if(nFoundppi>3) return -1;
2247  }else if(TMath::Abs(pdgdau)==2212){
2248  nProtons++;
2249  sumPxDau+=dau->Px();
2250  sumPyDau+=dau->Py();
2251  sumPzDau+=dau->Pz();
2252  arrayDauLab[nFoundppi++]=indDau;
2253  if(nFoundppi>3) return -1;
2254  }else if(TMath::Abs(pdgdau)==311 || TMath::Abs(pdgdau)==3122){
2255  codeV0=TMath::Abs(pdgdau);
2256  TParticle* v0=dau;
2257  AliMCParticle* mcV0=mcDau;
2258  if(codeV0==311){
2259  Int_t nK0Dau=mcDau->GetNDaughters();
2260  if(nK0Dau!=1) return -1;
2261  Int_t indK0s=mcDau->GetDaughterFirst();
2262  if(indK0s<0) return -1;
2263  v0=mcEvent->Particle(indK0s);
2264  mcV0=(AliMCParticle*)mcEvent->GetTrack(indK0s);
2265  if(!v0 || !mcV0) return -1;
2266  Int_t pdgK0sl=v0->GetPdgCode();
2267  codeV0=TMath::Abs(pdgK0sl);
2268  }
2269  Int_t nV0Dau=mcV0->GetNDaughters();
2270  if(nV0Dau!=2) return -1;
2271  Int_t indFirstV0Dau=mcV0->GetDaughterFirst();
2272  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
2273  Int_t indV0Dau=indFirstV0Dau+v0Dau;
2274  if(indV0Dau<0) return -1;
2275  TParticle* v0dau=mcEvent->Particle(indV0Dau);
2276  if(!v0dau) return -1;
2277  Int_t pdgv0dau=v0dau->GetPdgCode();
2278  if(TMath::Abs(pdgv0dau)==211){
2279  sumPxDau+=v0dau->Px();
2280  sumPyDau+=v0dau->Py();
2281  sumPzDau+=v0dau->Pz();
2282  nPions++;
2283  arrayDauLab[nFoundppi++]=indV0Dau;
2284  if(nFoundppi>3) return -1;
2285  }else if(TMath::Abs(pdgv0dau)==2212){
2286  sumPxDau+=v0dau->Px();
2287  sumPyDau+=v0dau->Py();
2288  sumPzDau+=v0dau->Pz();
2289  nProtons++;
2290  arrayDauLab[nFoundppi++]=indV0Dau;
2291  if(nFoundppi>3) return -1;
2292  }
2293  }
2294  }else{
2295  return -1;
2296  }
2297  }
2298  if(nPions!=2) return -1;
2299  if(nProtons!=1) return -1;
2300  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
2301  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
2302  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
2303  if(codeV0==310) return 1;
2304  else if(codeV0==3122) return 2;
2305  }
2306  return -1;
2307 
2308 }
2309 //____________________________________________________________________________
2310 Int_t AliVertexingHFUtils::CheckLcV0bachelorDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
2312 
2313  Int_t pdgD=mcPart->GetPdgCode();
2314  if(TMath::Abs(pdgD)!=4122) return -1;
2315 
2316  Int_t nDau=mcPart->GetNDaughters();
2317  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
2318  Int_t nPions=0;
2319  Int_t nProtons=0;
2320  Double_t sumPxDau=0.;
2321  Double_t sumPyDau=0.;
2322  Double_t sumPzDau=0.;
2323  Int_t nFoundppi=0;
2324 
2325  Int_t codeV0=-1;
2326  if(nDau==2){
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)==211){
2334  nPions++;
2335  sumPxDau+=dau->Px();
2336  sumPyDau+=dau->Py();
2337  sumPzDau+=dau->Pz();
2338  arrayDauLab[nFoundppi++]=indDau;
2339  if(nFoundppi>3) return -1;
2340  }else if(TMath::Abs(pdgdau)==2212){
2341  nProtons++;
2342  sumPxDau+=dau->Px();
2343  sumPyDau+=dau->Py();
2344  sumPzDau+=dau->Pz();
2345  arrayDauLab[nFoundppi++]=indDau;
2346  if(nFoundppi>3) return -1;
2347  }else if(TMath::Abs(pdgdau)==311 || TMath::Abs(pdgdau)==3122){
2348  codeV0=TMath::Abs(pdgdau);
2349  AliAODMCParticle* v0=dau;
2350  if(codeV0==311){
2351  Int_t nK0Dau=dau->GetNDaughters();
2352  if(nK0Dau!=1) return -1;
2353  Int_t indK0s=dau->GetDaughterLabel(0);
2354  if(indK0s<0) return -1;
2355  v0=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indK0s));
2356  if(!v0) return -1;
2357  Int_t pdgK0sl=v0->GetPdgCode();
2358  codeV0=TMath::Abs(pdgK0sl);
2359  }
2360  Int_t nV0Dau=v0->GetNDaughters();
2361  if(nV0Dau!=2) return -1;
2362  Int_t indFirstV0Dau=v0->GetDaughterLabel(0);
2363  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
2364  Int_t indV0Dau=indFirstV0Dau+v0Dau;
2365  if(indV0Dau<0) return -1;
2366  AliAODMCParticle* v0dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indV0Dau));
2367  if(!v0dau) return -1;
2368  Int_t pdgv0dau=v0dau->GetPdgCode();
2369  if(TMath::Abs(pdgv0dau)==211){
2370  sumPxDau+=v0dau->Px();
2371  sumPyDau+=v0dau->Py();
2372  sumPzDau+=v0dau->Pz();
2373  nPions++;
2374  arrayDauLab[nFoundppi++]=indV0Dau;
2375  if(nFoundppi>3) return -1;
2376  }else if(TMath::Abs(pdgv0dau)==2212){
2377  sumPxDau+=v0dau->Px();
2378  sumPyDau+=v0dau->Py();
2379  sumPzDau+=v0dau->Pz();
2380  nProtons++;
2381  arrayDauLab[nFoundppi++]=indV0Dau;
2382  if(nFoundppi>3) return -1;
2383  }
2384  }
2385  }else{
2386  return -1;
2387  }
2388  }
2389  if(nPions!=2) return -1;
2390  if(nProtons!=1) return -1;
2391  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
2392  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
2393  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
2394  if(codeV0==310) return 1;
2395  else if(codeV0==3122) return 2;
2396  }
2397  return -1;
2398 
2399 }
2400 
2401 //__________________________________xic______________________________________
2402 Int_t AliVertexingHFUtils::CheckXicXipipiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
2404 
2405  if(label<0) return -1;
2406  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
2407  TParticle* part = mcEvent->Particle(label);
2408  if(!part || !mcPart) return -1;
2409  Int_t pdgD=part->GetPdgCode();
2410  if(TMath::Abs(pdgD)!=4232) return -1;
2411 
2412  Int_t nDau=mcPart->GetNDaughters();
2413  if(nDau!=3) return -1;
2414 
2415  Int_t labelFirstDau = mcPart->GetDaughterFirst();
2416  Int_t nXi=0;
2417  Int_t nPions=0;
2418  Double_t sumPxDau=0.;
2419  Double_t sumPyDau=0.;
2420  Double_t sumPzDau=0.;
2421  Int_t nFoundXi=0;
2422 
2423  for(Int_t iDau=0; iDau<nDau; iDau++){
2424  Int_t indDau = labelFirstDau+iDau;
2425  if(indDau<0) return -1;
2426  TParticle* dau=mcEvent->Particle(indDau);
2427  if(!dau) return -1;
2428  Int_t pdgdau=dau->GetPdgCode();
2429  if(TMath::Abs(pdgdau)==3312){
2430  if(pdgD*pdgdau<0) return -1;
2431  sumPxDau+=dau->Px();
2432  sumPyDau+=dau->Py();
2433  sumPzDau+=dau->Pz();
2434  nXi++;
2435  arrayDauLab[nFoundXi++]=indDau;
2436 
2437  }
2438  if(TMath::Abs(pdgdau)==211){
2439  if(pdgD*pdgdau<0) return -1;
2440  nPions++;
2441  sumPxDau+=dau->Px();
2442  sumPyDau+=dau->Py();
2443  sumPzDau+=dau->Pz();
2444  arrayDauLab[nFoundXi++]=indDau;
2445  if(nFoundXi>3) return -1;
2446  }
2447  }
2448 
2449  if(nPions!=2) return -1;
2450  if(nXi!=1) return -1;
2451  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
2452  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
2453  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
2454  return 1;
2455 
2456 }
2457 
2458 //____________________________________________________________________________
2459 Int_t AliVertexingHFUtils::CheckBplusDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
2464 
2465  if(label<0) return -1;
2466  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
2467  TParticle* part = mcEvent->Particle(label);
2468  if(!part || !mcPart) return -1;
2469  Int_t pdgD=part->GetPdgCode();
2470  if(TMath::Abs(pdgD)!=521) return -1;
2471 
2472  Int_t nDau=mcPart->GetNDaughters();
2473  if(nDau!=2) return -1;
2474 
2475  Int_t labelFirstDau = mcPart->GetDaughterFirst();
2476  Int_t nKaons=0;
2477  Int_t nPions=0;
2478  Double_t sumPxDau=0.;
2479  Double_t sumPyDau=0.;
2480  Double_t sumPzDau=0.;
2481  Int_t nFoundKpi=0;
2482 
2483  for(Int_t iDau=0; iDau<nDau; iDau++){
2484  Int_t indDau = labelFirstDau+iDau;
2485  if(indDau<0) return -1;
2486  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
2487  TParticle* dau=mcEvent->Particle(indDau);
2488  if(!mcDau || !dau) return -1;
2489  Int_t pdgdau=dau->GetPdgCode();
2490  if(TMath::Abs(pdgdau)==421){
2491  Int_t nResDau=mcDau->GetNDaughters();
2492  if(nResDau!=2) return -1;
2493  Int_t indFirstResDau=mcDau->GetDaughterFirst();
2494  for(Int_t resDau=0; resDau<2; resDau++){
2495  Int_t indResDau=indFirstResDau+resDau;
2496  if(indResDau<0) return -1;
2497  TParticle* resdau=mcEvent->Particle(indResDau);
2498  if(!resdau) return -1;
2499  Int_t pdgresdau=resdau->GetPdgCode();
2500  if(TMath::Abs(pdgresdau)==321){
2501  if(pdgD*pdgresdau<0) return -1;
2502  sumPxDau+=resdau->Px();
2503  sumPyDau+=resdau->Py();
2504  sumPzDau+=resdau->Pz();
2505  nKaons++;
2506  arrayDauLab[nFoundKpi++]=indResDau;
2507  if(nFoundKpi>3) return -1;
2508  }
2509  if(TMath::Abs(pdgresdau)==211){
2510  if(pdgD*pdgresdau>0) return -1;
2511  sumPxDau+=resdau->Px();
2512  sumPyDau+=resdau->Py();
2513  sumPzDau+=resdau->Pz();
2514  nPions++;
2515  arrayDauLab[nFoundKpi++]=indResDau;
2516  if(nFoundKpi>3) return -1;
2517  }
2518  }
2519  }else if(TMath::Abs(pdgdau)==211){
2520  if(pdgD*pdgdau<0) return -1;
2521  nPions++;
2522  sumPxDau+=dau->Px();
2523  sumPyDau+=dau->Py();
2524  sumPzDau+=dau->Pz();
2525  arrayDauLab[nFoundKpi++]=indDau;
2526  if(nFoundKpi>3) return -1;
2527  }
2528  }
2529 
2530  if(nPions!=2) return -1;
2531  if(nKaons!=1) return -1;
2532  //Momentum conservation for several beauty decays not satisfied at gen. level in Upgrade MC's.
2533  //Fix implemented, loosening cut from 0.001 to 0.1. If >0.1, (-1*decay - 1) is returned.
2534  if(TMath::Abs(part->Px()-sumPxDau)>0.1) return -2;
2535  if(TMath::Abs(part->Py()-sumPyDau)>0.1) return -2;
2536  if(TMath::Abs(part->Pz()-sumPzDau)>0.1) return -2;
2537  return 1;
2538 
2539 }
2540 //____________________________________________________________________________
2541 Int_t AliVertexingHFUtils::CheckBplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
2546 
2547  Int_t pdgD=mcPart->GetPdgCode();
2548  if(TMath::Abs(pdgD)!=521) return -1;
2549 
2550  Int_t nDau=mcPart->GetNDaughters();
2551  if(nDau!=2) return -1;
2552 
2553  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
2554  Int_t nKaons=0;
2555  Int_t nPions=0;
2556  Double_t sumPxDau=0.;
2557  Double_t sumPyDau=0.;
2558  Double_t sumPzDau=0.;
2559  Int_t nFoundKpi=0;
2560 
2561  for(Int_t iDau=0; iDau<nDau; iDau++){
2562  Int_t indDau = labelFirstDau+iDau;
2563  if(indDau<0) return -1;
2564  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
2565  if(!dau) return -1;
2566  Int_t pdgdau=dau->GetPdgCode();
2567  if(TMath::Abs(pdgdau)==421){
2568  Int_t nResDau=dau->GetNDaughters();
2569  if(nResDau!=2) return -1;
2570  Int_t indFirstResDau=dau->GetDaughterLabel(0);
2571  for(Int_t resDau=0; resDau<2; resDau++){
2572  Int_t indResDau=indFirstResDau+resDau;
2573  if(indResDau<0) return -1;
2574  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
2575  if(!resdau) return -1;
2576  Int_t pdgresdau=resdau->GetPdgCode();
2577  if(TMath::Abs(pdgresdau)==321){
2578  if(pdgD*pdgresdau<0) return -1;
2579  sumPxDau+=resdau->Px();
2580  sumPyDau+=resdau->Py();
2581  sumPzDau+=resdau->Pz();
2582  nKaons++;
2583  arrayDauLab[nFoundKpi++]=indResDau;
2584  if(nFoundKpi>3) return -1;
2585  }
2586  if(TMath::Abs(pdgresdau)==211){
2587  if(pdgD*pdgresdau>0) return -1;
2588  sumPxDau+=resdau->Px();
2589  sumPyDau+=resdau->Py();
2590  sumPzDau+=resdau->Pz();
2591  nPions++;
2592  arrayDauLab[nFoundKpi++]=indResDau;
2593  if(nFoundKpi>3) return -1;
2594  }
2595  }
2596  }else if(TMath::Abs(pdgdau)==211){
2597  if(pdgD*pdgdau<0) return -1;
2598  nPions++;
2599  sumPxDau+=dau->Px();
2600  sumPyDau+=dau->Py();
2601  sumPzDau+=dau->Pz();
2602  arrayDauLab[nFoundKpi++]=indDau;
2603  if(nFoundKpi>3) return -1;
2604  }
2605  }
2606 
2607  if(nPions!=2) return -1;
2608  if(nKaons!=1) return -1;
2609  //Momentum conservation for several beauty decays not satisfied at gen. level in Upgrade MC's.
2610  //Fix implemented, loosening cut from 0.001 to 0.1. If >0.1, (-1*decay - 1) is returned.
2611  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.1) return -2;
2612  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.1) return -2;
2613  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.1) return -2;
2614  return 1;
2615 
2616 }
2617 //____________________________________________________________________________
2618 Int_t AliVertexingHFUtils::CheckB0toDminuspiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
2624 
2625  if(label<0) return -1;
2626  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
2627  TParticle* part = mcEvent->Particle(label);
2628  if(!part || !mcPart) return -1;
2629  Int_t pdgD=part->GetPdgCode();
2630  if(TMath::Abs(pdgD)!=511) return -1;
2631 
2632  Int_t nDau=mcPart->GetNDaughters();
2633  if(nDau!=2) return -1;
2634 
2635  Int_t labelFirstDau = mcPart->GetDaughterFirst();
2636  Int_t nKaons=0;
2637  Int_t nPions=0;
2638  Double_t sumPxDau=0.;
2639  Double_t sumPyDau=0.;
2640  Double_t sumPzDau=0.;
2641  Int_t nFoundKpi=0;
2642  Int_t decayB0 = -1;
2643 
2644  for(Int_t iDau=0; iDau<nDau; iDau++){
2645  Int_t indDau = labelFirstDau+iDau;
2646  if(indDau<0) return -1;
2647  TParticle* dau=mcEvent->Particle(indDau);
2648  if(!dau) return -1;
2649  Int_t pdgdau=dau->GetPdgCode();
2650  if(TMath::Abs(pdgdau)==411){
2652  Int_t labDauDplus[3] = {-1,-1,-1};
2653  Int_t decayDplus = CheckDplusDecay(mcEvent, indDau, labDauDplus);
2654 
2655  //Momentum conservation for several beauty decays not satisfied at gen. level in Upgrade MC's.
2656  //Fix implemented, to still select correct Dplus decay
2657  if(decayDplus==-2){
2658  AliMCParticle* mcDplus = (AliMCParticle*)mcEvent->GetTrack(indDau);
2659  Int_t labelFirstDauDplus = mcDplus->GetDaughterFirst();
2660  if(mcDplus->GetNDaughters() > 1){
2661  TParticle* dauDplus1 = mcEvent->Particle(labelFirstDauDplus);
2662  TParticle* dauDplus2 = mcEvent->Particle(labelFirstDauDplus+1);
2663  if(dauDplus1 && dauDplus2){
2664  Int_t pdgdauDplus1=dauDplus1->GetPdgCode();
2665  Int_t pdgdauDplus2=dauDplus2->GetPdgCode();
2666  if(TMath::Abs(pdgdauDplus1)==313 || TMath::Abs(pdgdauDplus2)==313) decayDplus=2;
2667  }
2668  }
2669  }
2670 
2671  if (decayDplus < 0 || labDauDplus[0] == -1) return -1;
2672  decayB0 = decayDplus;
2673  nPions+=2;
2674  nKaons++;
2675  for(Int_t iDplus = 0; iDplus < 3; iDplus++){
2676  TParticle* dauDplus=mcEvent->Particle(labDauDplus[iDplus]);
2677  sumPxDau+=dauDplus->Px();
2678  sumPyDau+=dauDplus->Py();
2679  sumPzDau+=dauDplus->Pz();
2680  arrayDauLab[nFoundKpi++]=labDauDplus[iDplus];
2681  }
2682  if(nFoundKpi>4) return -1;
2683  }else if(TMath::Abs(pdgdau)==211){
2684  //Temp fix for B0->D-+pi- and B0bar->D+pi- decays in ITS upgrade productions
2685  //if(pdgD*pdgdau>0) return -1;
2686  nPions++;
2687  sumPxDau+=dau->Px();
2688  sumPyDau+=dau->Py();
2689  sumPzDau+=dau->Pz();
2690  arrayDauLab[nFoundKpi++]=indDau;
2691  if(nFoundKpi>4) return -1;
2692  }
2693  }
2694 
2695  if(nPions!=3) return -1;
2696  if(nKaons!=1) return -1;
2697  //Momentum conservation for several beauty decays not satisfied at gen. level in Upgrade MC's.
2698  //Fix implemented, loosening cut from 0.001 to 0.1. If >0.1, (-1*decay - 1) is returned.
2699  Int_t decayB0temp = decayB0;
2700  if(TMath::Abs(part->Px()-sumPxDau)>0.1) decayB0temp = -1*decayB0 - 1;
2701  if(TMath::Abs(part->Py()-sumPyDau)>0.1) decayB0temp = -1*decayB0 - 1;
2702  if(TMath::Abs(part->Pz()-sumPzDau)>0.1) decayB0temp = -1*decayB0 - 1;
2703  decayB0 = decayB0temp;
2704  return decayB0;
2705 }
2706 //____________________________________________________________________________
2707 Int_t AliVertexingHFUtils::CheckB0toDminuspiDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
2713 
2714  Int_t pdgD=mcPart->GetPdgCode();
2715  if(TMath::Abs(pdgD)!=511) return -1;
2716 
2717  Int_t nDau=mcPart->GetNDaughters();
2718  if(nDau!=2) return -1;
2719 
2720  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
2721  Int_t nKaons=0;
2722  Int_t nPions=0;
2723  Double_t sumPxDau=0.;
2724  Double_t sumPyDau=0.;
2725  Double_t sumPzDau=0.;
2726  Int_t nFoundKpi=0;
2727  Int_t decayB0 = -1;
2728 
2729  for(Int_t iDau=0; iDau<nDau; iDau++){
2730  Int_t indDau = labelFirstDau+iDau;
2731  if(indDau<0) return -1;
2732  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
2733  if(!dau) return -1;
2734  Int_t pdgdau=dau->GetPdgCode();
2735  if(TMath::Abs(pdgdau)==411){
2737  Int_t labDauDplus[3] = {-1,-1,-1};
2738  Int_t decayDplus = CheckDplusDecay(arrayMC, dau, labDauDplus);
2739 
2740  //Momentum conservation for several beauty decays not satisfied at gen. level in Upgrade MC's.
2741  //Fix implemented, to still select correct Dplus decay
2742  if(decayDplus==-2){
2743  Int_t labelFirstDauDplus = dau->GetDaughterLabel(0);
2744  if(dau->GetNDaughters() > 1){
2745  AliAODMCParticle* dauDplus1=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labelFirstDauDplus));
2746  AliAODMCParticle* dauDplus2=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labelFirstDauDplus+1));
2747  if(dauDplus1 && dauDplus2){
2748  Int_t pdgdauDplus1=dauDplus1->GetPdgCode();
2749  Int_t pdgdauDplus2=dauDplus2->GetPdgCode();
2750  if(TMath::Abs(pdgdauDplus1)==313 || TMath::Abs(pdgdauDplus2)==313) decayDplus=2;
2751  }
2752  }
2753  }
2754 
2755  if (decayDplus < 0 || labDauDplus[0] == -1) return -1;
2756  decayB0 = decayDplus;
2757  nPions+=2;
2758  nKaons++;
2759  for(Int_t iDplus = 0; iDplus < 3; iDplus++){
2760  AliAODMCParticle* dauDplus=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDauDplus[iDplus]));
2761  sumPxDau+=dauDplus->Px();
2762  sumPyDau+=dauDplus->Py();
2763  sumPzDau+=dauDplus->Pz();
2764  arrayDauLab[nFoundKpi++]=labDauDplus[iDplus];
2765  }
2766  if(nFoundKpi>4) return -1;
2767  }else if(TMath::Abs(pdgdau)==211){
2768  //Temp fix for B0->D-+pi- and B0bar->D+pi- decays in ITS upgrade productions
2769  //if(pdgD*pdgdau>0) return -1;
2770  nPions++;
2771  sumPxDau+=dau->Px();
2772  sumPyDau+=dau->Py();
2773  sumPzDau+=dau->Pz();
2774  arrayDauLab[nFoundKpi++]=indDau;
2775  if(nFoundKpi>4) return -1;
2776  }
2777  }
2778 
2779  if(nPions!=3) return -1;
2780  if(nKaons!=1) return -1;
2781  //Momentum conservation for several beauty decays not satisfied at gen. level in Upgrade MC's.
2782  //Fix implemented, loosening cut from 0.001 to 0.1. If >0.1, (-1*decay - 1) is returned.
2783  Int_t decayB0temp = decayB0;
2784  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.1) decayB0temp = -1*decayB0 - 1;
2785  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.1) decayB0temp = -1*decayB0 - 1;
2786  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.1) decayB0temp = -1*decayB0 - 1;
2787  decayB0 = decayB0temp;
2788  return decayB0;
2789 }
2790 //____________________________________________________________________________
2791 Int_t AliVertexingHFUtils::CheckBsDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab, Bool_t ITS2UpgradeProd){
2797 
2798  if(label<0) return -1;
2799  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
2800  TParticle* part = mcEvent->Particle(label);
2801  if(!part || !mcPart) return -1;
2802  Int_t pdgD=part->GetPdgCode();
2803  if(TMath::Abs(pdgD)!=531) return -1;
2804 
2805  Int_t nDau=mcPart->GetNDaughters();
2806  if(nDau!=2) return -1;
2807 
2808  Int_t labelFirstDau = mcPart->GetDaughterFirst();
2809  Int_t nKaons=0;
2810  Int_t nPions=0;
2811  Double_t sumPxDau=0.;
2812  Double_t sumPyDau=0.;
2813  Double_t sumPzDau=0.;
2814  Int_t nFoundKpi=0;
2815  Int_t decayBs = -1;
2816 
2817  for(Int_t iDau=0; iDau<nDau; iDau++){
2818  Int_t indDau = labelFirstDau+iDau;
2819  if(indDau<0) return -1;
2820  TParticle* dau=mcEvent->Particle(indDau);
2821  if(!dau) return -1;
2822  Int_t pdgdau=dau->GetPdgCode();
2823  if(TMath::Abs(pdgdau)==431){
2824  //Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case, 4 for Ds->f0pi->KKpi
2825  Int_t labDauDs[3] = {-1,-1,-1};
2826  Int_t decayDs = CheckDsDecay(mcEvent, indDau, labDauDs);
2827 
2828  //Momentum conservation for several beauty decays not satisfied at gen. level in Upgrade MC's.
2829  //Fix implemented, to still select correct Ds decay
2830  if(decayDs==-2){
2831  AliMCParticle* mcDs = (AliMCParticle*)mcEvent->GetTrack(indDau);
2832  Int_t labelFirstDauDs = mcDs->GetDaughterFirst();
2833  if(mcDs->GetNDaughters() > 1){
2834  TParticle* dauDs1 = mcEvent->Particle(labelFirstDauDs);
2835  TParticle* dauDs2 = mcEvent->Particle(labelFirstDauDs+1);
2836  if(dauDs1 && dauDs2){
2837  Int_t pdgdauDs1=dauDs1->GetPdgCode();
2838  Int_t pdgdauDs2=dauDs2->GetPdgCode();
2839  if(TMath::Abs(pdgdauDs1)==313 || TMath::Abs(pdgdauDs2)==313) decayDs=2;
2840  else if(TMath::Abs(pdgdauDs1)==333 || TMath::Abs(pdgdauDs2)==333) decayDs=1;
2841  else if(TMath::Abs(pdgdauDs1)==9010221 || TMath::Abs(pdgdauDs2)==9010221) decayDs=4;
2842  else decayDs=3;
2843  }
2844  }
2845  }
2846 
2847  //In ITS2 Upgrade production, phi not "stored", so decay read as non-resonant
2848  if(decayDs==3 && ITS2UpgradeProd){
2849  TParticle* dauK1 = mcEvent->Particle(labDauDs[0]);
2850  TParticle* dauK2 = mcEvent->Particle(labDauDs[1]);
2851  TParticle* dauK3 = mcEvent->Particle(labDauDs[2]);
2852 
2853  TLorentzVector vK1, vK2, vKK;
2854  if(TMath::Abs(dauK1->GetPdgCode())==321 && TMath::Abs(dauK2->GetPdgCode())==321){
2855  dauK1->Momentum(vK1);
2856  dauK2->Momentum(vK2);
2857  } else if(TMath::Abs(dauK1->GetPdgCode())==321 && TMath::Abs(dauK3->GetPdgCode())==321){
2858  dauK1->Momentum(vK1);
2859  dauK3->Momentum(vK2);
2860  } else {
2861  dauK2->Momentum(vK1);
2862  dauK3->Momentum(vK2);
2863  }
2864  vKK = vK1 + vK2;
2865  //Small window around phi-mass, tag as Ds->phipi->KKpi if inside
2866  if(vKK.M() > 1.00 && vKK.M() < 1.04) decayDs = 1;
2867  }
2868 
2869  if (decayDs < 0 || labDauDs[0] == -1) return -1;
2870  decayBs = decayDs;
2871  nPions++;
2872  nKaons+=2;
2873  for(Int_t iDs = 0; iDs < 3; iDs++){
2874  TParticle* dauDs=mcEvent->Particle(labDauDs[iDs]);
2875  sumPxDau+=dauDs->Px();
2876  sumPyDau+=dauDs->Py();
2877  sumPzDau+=dauDs->Pz();
2878  arrayDauLab[nFoundKpi++]=labDauDs[iDs];
2879  }
2880  if(nFoundKpi>4) return -1;
2881  }else if(TMath::Abs(pdgdau)==211){
2882  //Temp fix for Bs->Ds+pi- and Bs->Ds-pi+ decays in ITS upgrade productions
2883  //if(pdgD*pdgdau>0) return -1;
2884  nPions++;
2885  sumPxDau+=dau->Px();
2886  sumPyDau+=dau->Py();
2887  sumPzDau+=dau->Pz();
2888  arrayDauLab[nFoundKpi++]=indDau;
2889  if(nFoundKpi>4) return -1;
2890  }
2891  }
2892 
2893  if(nPions!=2) return -1;
2894  if(nKaons!=2) return -1;
2895  //Momentum conservation for several beauty decays not satisfied at gen. level in Upgrade MC's.
2896  //Fix implemented, loosening cut from 0.001 to 0.1. If >0.1, (-1*decay - 1) is returned.
2897  Int_t decayBstemp = decayBs;
2898  if(TMath::Abs(part->Px()-sumPxDau)>0.1) decayBstemp = -1*decayBs - 1;
2899  if(TMath::Abs(part->Py()-sumPyDau)>0.1) decayBstemp = -1*decayBs - 1;
2900  if(TMath::Abs(part->Pz()-sumPzDau)>0.1) decayBstemp = -1*decayBs - 1;
2901  decayBs = decayBstemp;
2902  return decayBs;
2903 }
2904 //____________________________________________________________________________
2905 Int_t AliVertexingHFUtils::CheckBsDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab, Bool_t ITS2UpgradeProd){
2911 
2912  Int_t pdgD=mcPart->GetPdgCode();
2913  if(TMath::Abs(pdgD)!=531) return -1;
2914 
2915  Int_t nDau=mcPart->GetNDaughters();
2916  if(nDau!=2) return -1;
2917 
2918  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
2919  Int_t nKaons=0;
2920  Int_t nPions=0;
2921  Double_t sumPxDau=0.;
2922  Double_t sumPyDau=0.;
2923  Double_t sumPzDau=0.;
2924  Int_t nFoundKpi=0;
2925  Int_t decayBs = -1;
2926 
2927  for(Int_t iDau=0; iDau<nDau; iDau++){
2928  Int_t indDau = labelFirstDau+iDau;
2929  if(indDau<0) return -1;
2930  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
2931  if(!dau) return -1;
2932  Int_t pdgdau=dau->GetPdgCode();
2933  if(TMath::Abs(pdgdau)==431){
2934  //Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case, 4 for Ds->f0pi->KKpi
2935  Int_t labDauDs[3] = {-1,-1,-1};
2936  Int_t decayDs = CheckDsDecay(arrayMC, dau, labDauDs);
2937 
2938  //Momentum conservation for several beauty decays not satisfied at gen. level in Upgrade MC's.
2939  //Fix implemented, to still select correct Ds decay
2940  if(decayDs==-2){
2941  Int_t labelFirstDauDs = dau->GetDaughterLabel(0);
2942  if(dau->GetNDaughters() > 1){
2943  AliAODMCParticle* dauDs1=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labelFirstDauDs));
2944  AliAODMCParticle* dauDs2=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labelFirstDauDs+1));
2945  if(dauDs1 && dauDs2){
2946  Int_t pdgdauDs1=dauDs1->GetPdgCode();
2947  Int_t pdgdauDs2=dauDs2->GetPdgCode();
2948  if(TMath::Abs(pdgdauDs1)==313 || TMath::Abs(pdgdauDs2)==313) decayDs=2;
2949  else if(TMath::Abs(pdgdauDs1)==333 || TMath::Abs(pdgdauDs2)==333) decayDs=1;
2950  else if(TMath::Abs(pdgdauDs1)==9010221 || TMath::Abs(pdgdauDs2)==9010221) decayDs=4;
2951  else decayDs=3;
2952  }
2953  }
2954  }
2955 
2956  //In ITS2 Upgrade production, phi not "stored", so decay read as non-resonant
2957  if(decayDs==3 && ITS2UpgradeProd){
2958  AliAODMCParticle* dauK1=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDauDs[0]));
2959  AliAODMCParticle* dauK2=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDauDs[1]));
2960  AliAODMCParticle* dauK3=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDauDs[2]));
2961 
2962  TLorentzVector vK1, vK2, vKK;
2963  if(TMath::Abs(dauK1->GetPdgCode())==321 && TMath::Abs(dauK2->GetPdgCode())==321){
2964  dauK1->Momentum(vK1);
2965  dauK2->Momentum(vK2);
2966  } else if(TMath::Abs(dauK1->GetPdgCode())==321 && TMath::Abs(dauK3->GetPdgCode())==321){
2967  dauK1->Momentum(vK1);
2968  dauK3->Momentum(vK2);
2969  } else {
2970  dauK2->Momentum(vK1);
2971  dauK3->Momentum(vK2);
2972  }
2973  vKK = vK1 + vK2;
2974  //Small window around phi-mass, tag as Ds->phipi->KKpi if inside
2975  if(vKK.M() > 1.00 && vKK.M() < 1.04) decayDs = 1;
2976  }
2977 
2978  if (decayDs < 0 || labDauDs[0] == -1) return -1;
2979  decayBs = decayDs;
2980  nPions++;
2981  nKaons+=2;
2982  for(Int_t iDs = 0; iDs < 3; iDs++){
2983  AliAODMCParticle* dauDs=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDauDs[iDs]));
2984  sumPxDau+=dauDs->Px();
2985  sumPyDau+=dauDs->Py();
2986  sumPzDau+=dauDs->Pz();
2987  arrayDauLab[nFoundKpi++]=labDauDs[iDs];
2988  }
2989  if(nFoundKpi>4) return -1;
2990  }else if(TMath::Abs(pdgdau)==211){
2991  //Temp fix for Bs->Ds+pi- and Bs->Ds-pi+ decays in ITS upgrade productions
2992  //if(pdgD*pdgdau>0) return -1;
2993  nPions++;
2994  sumPxDau+=dau->Px();
2995  sumPyDau+=dau->Py();
2996  sumPzDau+=dau->Pz();
2997  arrayDauLab[nFoundKpi++]=indDau;
2998  if(nFoundKpi>4) return -1;
2999  }
3000  }
3001 
3002  if(nPions!=2) return -1;
3003  if(nKaons!=2) return -1;
3004  //Momentum conservation for several beauty decays not satisfied at gen. level in Upgrade MC's.
3005  //Fix implemented, loosening cut from 0.001 to 0.1. If >0.1, (-1*decay - 1) is returned.
3006  Int_t decayBstemp = decayBs;
3007  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.1) decayBstemp = -1*decayBs - 1;
3008  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.1) decayBstemp = -1*decayBs - 1;
3009  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.1) decayBstemp = -1*decayBs - 1;
3010  decayBs = decayBstemp;
3011  return decayBs;
3012 }
3013 //____________________________________________________________________________
3014 Int_t AliVertexingHFUtils::CheckLbDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
3020 
3021  if(label<0) return -1;
3022  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
3023  TParticle* part = mcEvent->Particle(label);
3024  if(!part || !mcPart) return -1;
3025  Int_t pdgD=part->GetPdgCode();
3026  if(TMath::Abs(pdgD)!=5122) return -1;
3027 
3028  Int_t nDau=mcPart->GetNDaughters();
3029  if(nDau!=2) return -1;
3030 
3031  Int_t labelFirstDau = mcPart->GetDaughterFirst();
3032  Int_t nKaons=0;
3033  Int_t nPions=0;
3034  Int_t nProtons=0;
3035  Double_t sumPxDau=0.;
3036  Double_t sumPyDau=0.;
3037  Double_t sumPzDau=0.;
3038  Int_t nFoundpKpi=0;
3039  Int_t decayLb = -1;
3040 
3041  for(Int_t iDau=0; iDau<nDau; iDau++){
3042  Int_t indDau = labelFirstDau+iDau;
3043  if(indDau<0) return -1;
3044  TParticle* dau=mcEvent->Particle(indDau);
3045  if(!dau) return -1;
3046  Int_t pdgdau=dau->GetPdgCode();
3047  if(TMath::Abs(pdgdau)==4122){
3048  Int_t labDauLc[3] = {-1,-1,-1};
3049  //Returns 1 for non-resonant decays and 2, 3 or 4 for resonant ones, -1 in other cases
3050  Int_t decayLc = CheckLcpKpiDecay(mcEvent, indDau, labDauLc);
3051 
3052  //Momentum conservation for several beauty decays not satisfied at gen. level in Upgrade MC's.
3053  //Fix implemented, to still select correct Lc decay
3054  if(decayLc==-2){
3055  AliMCParticle* mcLc = (AliMCParticle*)mcEvent->GetTrack(indDau);
3056  Int_t labelFirstDauLc = mcLc->GetDaughterFirst();
3057  if(mcLc->GetNDaughters() > 1){
3058  TParticle* dauLc1 = mcEvent->Particle(labelFirstDauLc);
3059  TParticle* dauLc2 = mcEvent->Particle(labelFirstDauLc+1);
3060  if(dauLc1 && dauLc2){
3061  Int_t pdgdauLc1=dauLc1->GetPdgCode();
3062  Int_t pdgdauLc2=dauLc2->GetPdgCode();
3063  if(TMath::Abs(pdgdauLc1)==313 || TMath::Abs(pdgdauLc2)==313) decayLc=2;
3064  else if(TMath::Abs(pdgdauLc1)==2224 || TMath::Abs(pdgdauLc2)==2224) decayLc=3;
3065  else if(TMath::Abs(pdgdauLc1)==3124 || TMath::Abs(pdgdauLc2)==3124) decayLc=4;
3066  else decayLc=1;
3067  }
3068  }
3069  }
3070 
3071  if (decayLc < 0 || labDauLc[0] == -1) return -1;
3072  decayLb = decayLc;
3073  nProtons++;
3074  nKaons++;
3075  nPions++;
3076  for(Int_t iLc = 0; iLc < 3; iLc++){
3077  TParticle* dauLc=mcEvent->Particle(labDauLc[iLc]);
3078  sumPxDau+=dauLc->Px();
3079  sumPyDau+=dauLc->Py();
3080  sumPzDau+=dauLc->Pz();
3081  arrayDauLab[nFoundpKpi++]=labDauLc[iLc];
3082  }
3083  if(nFoundpKpi>4) return -1;
3084  }else if(TMath::Abs(pdgdau)==211){
3085  //Temp fix for Lb->Lc+pi- and Lb->Lc-pi+ decays in ITS upgrade productions
3086  //if(pdgD*pdgdau>0) return -1;
3087  nPions++;
3088  sumPxDau+=dau->Px();
3089  sumPyDau+=dau->Py();
3090  sumPzDau+=dau->Pz();
3091  arrayDauLab[nFoundpKpi++]=indDau;
3092  if(nFoundpKpi>4) return -1;
3093  }
3094  }
3095 
3096  if(nProtons!=1) return -1;
3097  if(nKaons!=1) return -1;
3098  if(nPions!=2) return -1;
3099  //Momentum conservation for several beauty decays not satisfied at gen. level in Upgrade MC's.
3100  //Fix implemented, loosening cut from 0.001 to 0.1. If >0.1, (-1*decay - 1) is returned.
3101  Int_t decayLbtemp = decayLb;
3102  if(TMath::Abs(part->Px()-sumPxDau)>0.1) decayLbtemp = -1*decayLb - 1;
3103  if(TMath::Abs(part->Py()-sumPyDau)>0.1) decayLbtemp = -1*decayLb - 1;
3104  if(TMath::Abs(part->Pz()-sumPzDau)>0.1) decayLbtemp = -1*decayLb - 1;
3105  decayLb = decayLbtemp;
3106  return decayLb;
3107 }
3108 //____________________________________________________________________________
3109 Int_t AliVertexingHFUtils::CheckLbDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
3115 
3116  Int_t pdgD=mcPart->GetPdgCode();
3117  if(TMath::Abs(pdgD)!=5122) return -1;
3118 
3119  Int_t nDau=mcPart->GetNDaughters();
3120  if(nDau!=2) return -1;
3121 
3122  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
3123  Int_t nKaons=0;
3124  Int_t nPions=0;
3125  Int_t nProtons=0;
3126  Double_t sumPxDau=0.;
3127  Double_t sumPyDau=0.;
3128  Double_t sumPzDau=0.;
3129  Int_t nFoundpKpi=0;
3130  Int_t decayLb = -1;
3131 
3132  for(Int_t iDau=0; iDau<nDau; iDau++){
3133  Int_t indDau = labelFirstDau+iDau;
3134  if(indDau<0) return -1;
3135  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
3136  if(!dau) return -1;
3137  Int_t pdgdau=dau->GetPdgCode();
3138  if(TMath::Abs(pdgdau)==4122){
3139  Int_t labDauLc[3] = {-1,-1,-1};
3140  //Returns 1 for non-resonant decays and 2, 3 or 4 for resonant ones, -1 in other cases
3141  Int_t decayLc = CheckLcpKpiDecay(arrayMC, dau, labDauLc);
3142 
3143  //Momentum conservation for several beauty decays not satisfied at gen. level in Upgrade MC's.
3144  //Fix implemented, to still select correct Lc decay
3145  if(decayLc==-2){
3146  Int_t labelFirstDauLc = dau->GetDaughterLabel(0);
3147  if(dau->GetNDaughters() > 1){
3148  AliAODMCParticle* dauLc1=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labelFirstDauLc));
3149  AliAODMCParticle* dauLc2=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labelFirstDauLc+1));
3150  if(dauLc1 && dauLc2){
3151  Int_t pdgdauLc1=dauLc1->GetPdgCode();
3152  Int_t pdgdauLc2=dauLc2->GetPdgCode();
3153  if(TMath::Abs(pdgdauLc1)==313 || TMath::Abs(pdgdauLc2)==313) decayLc=2;
3154  else if(TMath::Abs(pdgdauLc1)==2224 || TMath::Abs(pdgdauLc2)==2224) decayLc=3;
3155  else if(TMath::Abs(pdgdauLc1)==3124 || TMath::Abs(pdgdauLc2)==3124) decayLc=4;
3156  else decayLc=1;
3157  }
3158  }
3159  }
3160 
3161  if (decayLc < 0 || labDauLc[0] == -1) return -1;
3162  decayLb = decayLc;
3163  nProtons++;
3164  nKaons++;
3165  nPions++;
3166  for(Int_t iLc = 0; iLc < 3; iLc++){
3167  AliAODMCParticle* dauLc=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDauLc[iLc]));
3168  sumPxDau+=dauLc->Px();
3169  sumPyDau+=dauLc->Py();
3170  sumPzDau+=dauLc->Pz();
3171  arrayDauLab[nFoundpKpi++]=labDauLc[iLc];
3172  }
3173  if(nFoundpKpi>4) return -1;
3174  }else if(TMath::Abs(pdgdau)==211){
3175  //Temp fix for Lb->Lc+pi- and Lb->Lc-pi+ decays in ITS upgrade productions
3176  //if(pdgD*pdgdau>0) return -1;
3177  nPions++;
3178  sumPxDau+=dau->Px();
3179  sumPyDau+=dau->Py();
3180  sumPzDau+=dau->Pz();
3181  arrayDauLab[nFoundpKpi++]=indDau;
3182  if(nFoundpKpi>4) return -1;
3183  }
3184  }
3185 
3186  if(nProtons!=1) return -1;
3187  if(nKaons!=1) return -1;
3188  if(nPions!=2) return -1;
3189  //Momentum conservation for several beauty decays not satisfied at gen. level in Upgrade MC's.
3190  //Fix implemented, loosening cut from 0.001 to 0.1. If >0.1, (-1*decay - 1) is returned.
3191  Int_t decayLbtemp = decayLb;
3192  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.1) decayLbtemp = -1*decayLb - 1;
3193  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.1) decayLbtemp = -1*decayLb - 1;
3194  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.1) decayLbtemp = -1*decayLb - 1;
3195  decayLb = decayLbtemp;
3196  return decayLb;
3197 }
3198 //________________________________________________________________________
3200  Double_t etaMin, Double_t etaMax,
3202  Int_t filtbit1, Int_t filtbit2,
3203  Int_t minMult){
3205 
3206  Int_t nTracks=aod->GetNumberOfTracks();
3207  Int_t nSelTracks=0;
3208 
3209  Double_t sumpt=0.;
3210  Double_t s00=0.;
3211  Double_t s01=0.;
3212  Double_t s11=0.;
3213  if(ptMin<0.) ptMin=0.;
3214 
3215  for(Int_t it=0; it<nTracks; it++) {
3216  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
3217  if(!tr) continue;
3218  Float_t eta = tr->Eta();
3219  Float_t pt = tr->Pt();
3220  Float_t phi = tr->Phi();
3221  if(eta<etaMin || eta>etaMax) continue;
3222  if(pt<ptMin || pt>ptMax) continue;
3223  Bool_t fb1 = tr->TestFilterBit(filtbit1);
3224  Bool_t fb2 = tr->TestFilterBit(filtbit2);
3225  Bool_t tpcRefit=tr->GetStatus() & AliAODTrack::kTPCrefit;
3226  if(filtbit1==1 && !tpcRefit) fb1=kFALSE;
3227  if(filtbit2==1 && !tpcRefit) fb2=kFALSE;
3228  if( !(fb1 || fb2) ) continue;
3229  Double_t px=pt*TMath::Cos(phi);
3230  Double_t py=pt*TMath::Sin(phi);
3231  s00 += (px * px)/pt;
3232  s01 += (py * px)/pt;
3233  s11 += (py * py)/pt;
3234  nSelTracks++;
3235  sumpt+=pt;
3236  }
3237 
3238  if(nSelTracks<minMult) return -0.5;
3239 
3240  if(sumpt>0.){
3241  s00/=sumpt;
3242  s01/=sumpt;
3243  s11/=sumpt;
3244  }else return -0.5;
3245 
3246  Double_t sphericity = -10;
3247  Double_t lambda1=((s00+s11)+TMath::Sqrt((s00+s11)*(s00+s11)-4*(s00*s11-s01*s01)))/2.;
3248  Double_t lambda2=((s00+s11)-TMath::Sqrt((s00+s11)*(s00+s11)-4*(s00*s11-s01*s01)))/2.;
3249  if(TMath::Abs(lambda2)<0.00001 && TMath::Abs(lambda1)<0.00001) sphericity=0;
3250  if(TMath::Abs(lambda1+lambda2)>0.000001) sphericity=2*TMath::Min(lambda1,lambda2)/(lambda1+lambda2);
3251  return sphericity;
3252 
3253 }
3254 
3255 //________________________________________________________________________
3257  Double_t &spherocity, Double_t &phiRef,
3258  Double_t etaMin, Double_t etaMax,
3260  Int_t filtbit1, Int_t filtbit2,
3261  Int_t minMult, Double_t phiStepSizeDeg,
3262  Int_t nTrksToSkip, Int_t* idToSkip
3263  ){
3265 
3266  Int_t nTracks=aod->GetNumberOfTracks();
3267  Int_t nSelTracks=0;
3268 
3269  Double_t* ptArr=new Double_t[nTracks];
3270  Double_t* phiArr=new Double_t[nTracks];
3271  Double_t sumpt=0.;
3272 
3273  for(Int_t it=0; it<nTracks; it++) {
3274  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
3275  if(!tr) continue;
3276  Float_t eta = tr->Eta();
3277  Float_t pt = tr->Pt();
3278  Float_t phi = tr->Phi();
3279  if(eta<etaMin || eta>etaMax) continue;
3280  if(pt<ptMin || pt>ptMax) continue;
3281  if(nTrksToSkip>0 && idToSkip){
3282  Int_t trid = (Int_t)tr->GetID();
3283  Bool_t keep=kTRUE;
3284  for(Int_t jt=0; jt<nTrksToSkip; jt++){
3285  if(trid==idToSkip[jt]) keep=kFALSE;
3286  }
3287  if(!keep) continue;
3288  }
3289  Bool_t fb1 = tr->TestFilterBit(filtbit1);
3290  Bool_t fb2 = tr->TestFilterBit(filtbit2);
3291  Bool_t tpcRefit=tr->GetStatus() & AliAODTrack::kTPCrefit;
3292  if(filtbit1==1 && !tpcRefit) fb1=kFALSE;
3293  if(filtbit2==1 && !tpcRefit) fb2=kFALSE;
3294  if( !(fb1 || fb2) ) continue;
3295  ptArr[nSelTracks]=pt;
3296  phiArr[nSelTracks]=phi;
3297  nSelTracks++;
3298  sumpt+=pt;
3299  }
3300 
3301  if(nSelTracks<minMult){spherocity = -0.5; return;}
3302 
3303  //Getting thrust
3304  spherocity=2.;
3305  for(Int_t i=0; i<360/phiStepSizeDeg; ++i){
3306  Double_t phistep=TMath::Pi()*(Double_t)i*phiStepSizeDeg/180.;
3307  Double_t nx=TMath::Cos(phistep);
3308  Double_t ny=TMath::Sin(phistep);
3309  Double_t numer=0.;
3310  for(Int_t j=0; j<nSelTracks; ++j){
3311  Double_t pxA=ptArr[j]*TMath::Cos(phiArr[j]); // x component of an unitary vector n
3312  Double_t pyA=ptArr[j]*TMath::Sin(phiArr[j]); // y component of an unitary vector n
3313  numer+=TMath::Abs(ny*pxA - nx*pyA);
3314  }
3315  Double_t pFull=numer*numer/(sumpt*sumpt);
3316  if(pFull<spherocity){
3317  spherocity=pFull; // minimization;
3318  phiRef=phistep;
3319  }
3320  }
3321 
3322  delete [] ptArr;
3323  delete [] phiArr;
3324 
3325  spherocity*=(TMath::Pi()*TMath::Pi()/4.);
3326  return;
3327 
3328 }
3329 //________________________________________________________________________
3331  Double_t &spherocity, Double_t &phiRef,
3332  Double_t etaMin, Double_t etaMax,
3334  Int_t minMult, Double_t phiStepSizeDeg){
3335 
3337 
3338  Int_t nParticles=arrayMC->GetEntriesFast();
3339  Int_t nSelParticles=0;
3340 
3341  Double_t* ptArr=new Double_t[nParticles];
3342  Double_t* phiArr=new Double_t[nParticles];
3343  Double_t sumpt=0.;
3344 
3345  for(Int_t ip=0; ip<nParticles; ip++) {
3346  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(ip);
3347  if(!part) continue;
3348  Float_t eta = part->Eta();
3349  Float_t pt = part->Pt();
3350  Float_t phi = part->Phi();
3351  Int_t charge = part->Charge();
3352  Bool_t isPhysPrim = part->IsPhysicalPrimary();
3353  if(!isPhysPrim) continue;
3354  if(charge==0) continue;
3355  if(eta<etaMin || eta>etaMax) continue;
3356  if(pt<ptMin || pt>ptMax) continue;
3357 
3358  ptArr[nSelParticles]=pt;
3359  phiArr[nSelParticles]=phi;
3360  nSelParticles++;
3361  sumpt+=pt;
3362  }
3363 
3364  if(nSelParticles<minMult){spherocity = -0.5; return;}
3365 
3366  //Getting thrust
3367  spherocity=2.;
3368  for(Int_t i=0; i<360/phiStepSizeDeg; ++i){
3369  Double_t phistep=TMath::Pi()*(Double_t)i*phiStepSizeDeg/180.;
3370  Double_t nx=TMath::Cos(phistep);
3371  Double_t ny=TMath::Sin(phistep);
3372  Double_t numer=0.;
3373  for(Int_t j=0; j<nSelParticles; ++j){
3374  Double_t pxA=ptArr[j]*TMath::Cos(phiArr[j]); // x component of an unitary vector n
3375  Double_t pyA=ptArr[j]*TMath::Sin(phiArr[j]); // y component of an unitary vector n
3376  numer+=TMath::Abs(ny*pxA - nx*pyA);
3377  }
3378  Double_t pFull=numer*numer/(sumpt*sumpt);
3379  if(pFull<spherocity){
3380  spherocity=pFull; // minimization;
3381  phiRef=phistep;
3382  }
3383  }
3384 
3385  delete [] ptArr;
3386  delete [] phiArr;
3387 
3388  spherocity*=(TMath::Pi()*TMath::Pi()/4.);
3389  return;
3390 
3391 }
3392 
3393 //________________________________________________________________________
3400 
3401  Int_t nBinOrig=hOrig->GetNbinsX();
3402  Int_t firstBinOrig=1;
3403  Int_t lastBinOrig=nBinOrig;
3404  Int_t nBinOrigUsed=nBinOrig;
3405  Int_t nBinFinal=nBinOrig/reb;
3406  if(firstUse>=1){
3407  firstBinOrig=firstUse;
3408  nBinFinal=(nBinOrig-firstUse+1)/reb;
3409  nBinOrigUsed=nBinFinal*reb;
3410  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
3411  }else{
3412  Int_t exc=nBinOrigUsed%reb;
3413  if(exc!=0){
3414  nBinOrigUsed-=exc;
3415  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
3416  }
3417  }
3418 
3419  printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);
3420  Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);
3421  Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);
3422  TH1D* hRebin=new TH1D(Form("%s-rebin",hOrig->GetName()),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);
3423  Int_t lastSummed=firstBinOrig-1;
3424  for(Int_t iBin=1;iBin<=nBinFinal; iBin++){
3425  Float_t sum=0.;
3426  Float_t sume2=0.;
3427  for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){
3428  sum+=hOrig->GetBinContent(lastSummed+1);
3429  sume2+=(hOrig->GetBinError(lastSummed+1)*hOrig->GetBinError(lastSummed+1));
3430  lastSummed++;
3431  }
3432  hRebin->SetBinContent(iBin,sum);
3433  hRebin->SetBinError(iBin,TMath::Sqrt(sume2));
3434  }
3435  return hRebin;
3436 }
3437 //________________________________________________________________________
3440 
3441  Int_t binmin=TMath::Max(1,hData->FindBin(hMC->GetXaxis()->GetXmin()));
3442  Bool_t found=kFALSE;
3443  Int_t binminD=-1;
3444  Int_t binminMC=-1;
3445  for(Int_t j=binmin; j<hData->GetNbinsX(); j++){
3446  if(found) break;
3447  for(Int_t k=1; k<hMC->GetNbinsX(); k++){
3448  Double_t delta=TMath::Abs(hMC->GetBinLowEdge(k)-hData->GetBinLowEdge(j));
3449  if(delta<0.0001){
3450  found=kTRUE;
3451  binminMC=k;
3452  binminD=j;
3453  }
3454  if(found) break;
3455  }
3456  }
3457  Int_t binmax=TMath::Min(hData->GetNbinsX(),hData->FindBin(hMC->GetXaxis()->GetXmax()*0.99999));
3458  found=kFALSE;
3459  Int_t binmaxD=-1;
3460  Int_t binmaxMC=-1;
3461  for(Int_t j=binmax; j>1; j--){
3462  if(found) break;
3463  for(Int_t k=hMC->GetNbinsX(); k>400; k--){
3464  Double_t delta=TMath::Abs(hMC->GetBinLowEdge(k+1)-hData->GetBinLowEdge(j+1));
3465  if(delta<0.0001){
3466  found=kTRUE;
3467  binmaxMC=k;
3468  binmaxD=j;
3469  }
3470  if(found) break;
3471  }
3472  }
3473 
3474  Double_t min=hData->GetBinLowEdge(binminD);
3475  Double_t max=hData->GetBinLowEdge(binmaxD)+hData->GetBinWidth(binmaxD);
3476  Double_t minMC=hMC->GetBinLowEdge(binminMC);
3477  Double_t maxMC=hMC->GetBinLowEdge(binmaxMC)+hMC->GetBinWidth(binmaxMC);
3478  Double_t width=hData->GetBinWidth(binminD);
3479  Double_t widthMC=hMC->GetBinWidth(binminMC);
3480 
3481  if(TMath::Abs(minMC-min)>0.0001*min || TMath::Abs(maxMC-max)>0.0001*max){
3482  printf("Cannot adapt range and rebin histo:\n");
3483  printf("Range for data histo: %f-%f GeV/c2 bins %d-%d width=%f\n",min,max,binminD,binmaxD,width);
3484  printf("Range for reflection histo: %f-%f GeV/c2 bins %d-%d width=%f\n",minMC,maxMC,binminMC,binmaxMC,widthMC);
3485  return 0x0;
3486  }
3487 
3488  Double_t rebin=width/widthMC;
3489  if(TMath::Abs(rebin-TMath::Nint(rebin))>0.001){
3490  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)));
3491  return 0x0;
3492  }
3493 
3494  Int_t nBinsNew=binmaxD-binminD+1;
3495  TH1 *hOut;
3496  TString stype=hMC->ClassName();
3497  if(stype.Contains("TH1F")){
3498  hOut=new TH1F(Form("%s-rebinned",hMC->GetName()),hMC->GetTitle(),nBinsNew,min,max);
3499  }else if(stype.Contains("TH1D")){
3500  hOut=new TH1D(Form("%s-rebinned",hMC->GetName()),hMC->GetTitle(),nBinsNew,min,max);
3501  }else{
3502  printf("Wrong type %s\n",stype.Data());
3503  return 0x0;
3504  }
3505 
3506  for(Int_t j=1; j<=hMC->GetNbinsX(); j++){
3507  Double_t m=hMC->GetBinCenter(j);
3508  Int_t binFin=hOut->FindBin(m);
3509  if(binFin>=1 && binFin<=nBinsNew){
3510  hOut->AddBinContent(binFin,hMC->GetBinContent(j));
3511  }
3512  }
3513  return hOut;
3514 }
3515 
3516 //___________________________________________________________________________________//
3517 //method that performs simultaneus fit of in-plane and out-of-plane inv-mass spectra
3518 ROOT::Fit::FitResult AliVertexingHFUtils::DoInPlaneOutOfPlaneSimultaneusFit(AliHFInvMassFitter *&massfitterInPlane, AliHFInvMassFitter *&massfitterOutOfPlane,
3519  TH1F* hMassInPlane, TH1F* hMassOutOfPlane, Double_t MinMass, Double_t MaxMass,
3520  Double_t massD, vector<UInt_t> commonpars) {
3521 
3522  cout << "\nIn-plane - out-of-plane simultaneus fit" << endl;
3523  cout << "\nIndependent prefits" << endl;
3524 
3525  //prefits to initialise parameters
3526  massfitterInPlane->SetUseLikelihoodFit();
3527  massfitterInPlane->SetInitialGaussianMean(massD);
3528  massfitterInPlane->SetInitialGaussianSigma(0.010);
3529  massfitterInPlane->MassFitter(kFALSE);
3530 
3531  massfitterOutOfPlane->SetUseLikelihoodFit();
3532  massfitterOutOfPlane->SetInitialGaussianMean(massD);
3533  massfitterOutOfPlane->SetInitialGaussianSigma(0.010);
3534  massfitterOutOfPlane->MassFitter(kFALSE);
3535 
3536  ROOT::Math::WrappedMultiTF1 wfInPlane(*(massfitterInPlane->GetMassFunc()),1);
3537  ROOT::Math::WrappedMultiTF1 wfOutOfPlane(*(massfitterOutOfPlane->GetMassFunc()),1);
3538 
3539  // set data options and ranges
3540  ROOT::Fit::DataOptions opt;
3541  ROOT::Fit::DataRange rangeMass; //same range for two functions
3542 
3543  rangeMass.SetRange(MinMass,MaxMass);
3544  ROOT::Fit::BinData dataInPlane(opt,rangeMass);
3545  ROOT::Fit::FillData(dataInPlane, hMassInPlane);
3546  ROOT::Fit::BinData dataOutOfPlane(opt,rangeMass);
3547  ROOT::Fit::FillData(dataOutOfPlane, hMassOutOfPlane);
3548 
3549  //define the 2 chi squares
3550  ROOT::Fit::Chi2Function chi2InPlane(dataInPlane, wfInPlane);
3551  ROOT::Fit::Chi2Function chi2OutOfPlane(dataOutOfPlane, wfOutOfPlane);
3552 
3553  //define the global chi square and get initial parameters from prefits
3554  const Int_t npars = massfitterInPlane->GetMassFunc()->GetNpar();
3555  const UInt_t ncommonpars = commonpars.size();
3556  Int_t nparsBkg = massfitterInPlane->GetBackgroundRecalcFunc()->GetNpar();
3557  Int_t nparsSgn = massfitterInPlane->GetSignalFunc()->GetNpar();
3558  Int_t nparsSecPeak = 0, nparsRefl = 0;
3559  if(massfitterInPlane->GetSecondPeakFunc())
3560  nparsSecPeak = massfitterInPlane->GetSecondPeakFunc()->GetNpar();
3561  if(massfitterInPlane->GetReflFunc())
3562  nparsRefl = massfitterInPlane->GetReflFunc()->GetNpar();
3563 
3564  GlobalInOutOfPlaneChi2 globalChi2(chi2InPlane, chi2OutOfPlane, npars, commonpars);
3565 
3566  vector<UInt_t>::iterator iter;
3567  vector<Double_t> initpars;
3568  for(Int_t iPar=0; iPar<2*npars; iPar++) {
3569  if(iPar<npars) { //in-plane
3570  initpars.push_back(massfitterInPlane->GetMassFunc()->GetParameter(iPar));
3571  }
3572  else { //out-of-plane
3573  iter = find(commonpars.begin(),commonpars.end(),iPar-npars);
3574  if(iter!=commonpars.end()) continue;
3575  else {
3576  initpars.push_back(massfitterOutOfPlane->GetMassFunc()->GetParameter(iPar-npars));
3577  }
3578  }
3579  }
3580 
3581  //define fitter and fit
3582  ROOT::Fit::Fitter simulfitter;
3583  simulfitter.Config().SetParamsSettings(npars*2-ncommonpars,initpars.data()); //set initial parameters from prefits
3584  if(nparsRefl>0) { //fix S/R
3585  simulfitter.Config().ParSettings(nparsBkg+nparsSgn+nparsSecPeak).Fix();
3586 
3587  iter = find(commonpars.begin(),commonpars.end(),npars+nparsBkg+nparsSgn+nparsSecPeak);
3588  if(iter==commonpars.end()) //if not included in common pars, need to be fixed also for out-of-plane func
3589  simulfitter.Config().ParSettings(npars+nparsBkg+nparsSgn+nparsSecPeak-ncommonpars).Fix();
3590  }
3591 
3592  simulfitter.Config().MinimizerOptions().SetPrintLevel(0);
3593  simulfitter.Config().SetMinimizer("Minuit2","Migrad");
3594  simulfitter.FitFCN(npars*2-ncommonpars,globalChi2,0,dataInPlane.Size()+dataOutOfPlane.Size(),kFALSE);
3595  ROOT::Fit::FitResult result = simulfitter.Result();
3596  cout << "\nSimultaneus fit" << endl;
3597  result.Print(cout);
3598 
3599  //Set new parameters to functions of mass fitters
3600  Int_t ncommonparsused = 0;
3601  for(Int_t iPar=0; iPar<npars; iPar++) {
3602  massfitterInPlane->GetMassFunc()->SetParameter(iPar,result.Parameter(iPar));
3603  massfitterInPlane->GetMassFunc()->SetParError(iPar,result.ParError(iPar));
3604  iter = find(commonpars.begin(),commonpars.end(),iPar);
3605  if(iter!=commonpars.end()) { //is common parameter
3606  massfitterOutOfPlane->GetMassFunc()->SetParameter(iPar,result.Parameter(iPar));
3607  massfitterOutOfPlane->GetMassFunc()->SetParError(iPar,result.ParError(iPar));
3608  ncommonparsused++;
3609  }
3610  else {
3611  massfitterOutOfPlane->GetMassFunc()->SetParameter(iPar,result.Parameter(iPar+npars-ncommonparsused));
3612  massfitterOutOfPlane->GetMassFunc()->SetParError(iPar,result.ParError(iPar+npars-ncommonparsused));
3613  }
3614 
3615  if(iPar < nparsBkg) { //bkg
3616  massfitterInPlane->GetBackgroundRecalcFunc()->SetParameter(iPar,massfitterInPlane->GetMassFunc()->GetParameter(iPar));
3617  massfitterInPlane->GetBackgroundRecalcFunc()->SetParError(iPar,massfitterInPlane->GetMassFunc()->GetParError(iPar));
3618  massfitterOutOfPlane->GetBackgroundRecalcFunc()->SetParameter(iPar,massfitterOutOfPlane->GetMassFunc()->GetParameter(iPar));
3619  massfitterOutOfPlane->GetBackgroundRecalcFunc()->SetParError(iPar,massfitterOutOfPlane->GetMassFunc()->GetParError(iPar));
3620  }
3621  else if(iPar >= nparsBkg && iPar < nparsBkg+nparsSgn){ //signal
3622  massfitterInPlane->GetSignalFunc()->SetParameter(iPar-nparsBkg,massfitterInPlane->GetMassFunc()->GetParameter(iPar));
3623  massfitterInPlane->GetSignalFunc()->SetParError(iPar-nparsBkg,massfitterInPlane->GetMassFunc()->GetParError(iPar));
3624  massfitterOutOfPlane->GetSignalFunc()->SetParameter(iPar-nparsBkg,massfitterOutOfPlane->GetMassFunc()->GetParameter(iPar));
3625  massfitterOutOfPlane->GetSignalFunc()->SetParError(iPar-nparsBkg,massfitterOutOfPlane->GetMassFunc()->GetParError(iPar));
3626  }
3627  else if(iPar >= nparsBkg && iPar < nparsBkg+nparsSgn){ //signal
3628  massfitterInPlane->GetSignalFunc()->SetParameter(iPar-nparsBkg,massfitterInPlane->GetMassFunc()->GetParameter(iPar));
3629  massfitterInPlane->GetSignalFunc()->SetParError(iPar-nparsBkg,massfitterInPlane->GetMassFunc()->GetParError(iPar));
3630  massfitterOutOfPlane->GetSignalFunc()->SetParameter(iPar-nparsBkg,massfitterOutOfPlane->GetMassFunc()->GetParameter(iPar));
3631  massfitterOutOfPlane->GetSignalFunc()->SetParError(iPar-nparsBkg,massfitterOutOfPlane->GetMassFunc()->GetParError(iPar));
3632  }
3633  }
3634 
3635  cout << "\n" << endl;
3636  return result;
3637 }
3638 
3639 //________________________________________________________________________
3641  // Compute maximum difference between observed and expected impact parameter of the candidate prongs
3642  Double_t dd0max = 0;
3643  UInt_t fNProngsCand = static_cast<UInt_t>(cand->GetNProngs());
3644  for (UInt_t iProng = 0; iProng < fNProngsCand; iProng++) {
3645  Double_t d0diff, errd0diff;
3646  cand->Getd0MeasMinusExpProng(iProng, bfield, d0diff, errd0diff);
3647  Double_t normdd0 = d0diff / errd0diff;
3648  if (iProng == 0 || TMath::Abs(normdd0) > TMath::Abs(dd0max))
3649  dd0max = normdd0;
3650  }
3651  return dd0max;
3652 }
3653 
3654 //______________________________________________________________________
3656  // Combine TPC and TOF nsigma (sum in quadrature + special cases)
3657  if (nsigmaTPC > -998. && nsigmaTOF > -998.)
3658  return TMath::Sqrt((nsigmaTPC * nsigmaTPC + nsigmaTOF * nsigmaTOF) / 2);
3659  else if (nsigmaTPC > -998. && nsigmaTOF < -998.)
3660  return TMath::Abs(nsigmaTPC);
3661  else if (nsigmaTPC < -998. && nsigmaTOF > -998.)
3662  return TMath::Abs(nsigmaTOF);
3663  else
3664  return -999.;
3665 }
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
void Getd0MeasMinusExpProng(Int_t ip, Double_t magf, Double_t &d0diff, Double_t &errd0diff) const
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
Int_t MassFitter(Bool_t draw=kTRUE)
static Double_t CombineNsigmaTPCTOF(Double_t nsigmaTPC, Double_t nsigmaTOF)
static TString GetGenerator(Int_t label, AliAODMCHeader *header)
AliHFInvMassFitter class for the fit of invariant mass distribution of charm hadrons.
static Bool_t IsTrackFromBeauty(AliAODTrack *tr, TClonesArray *arrayMC)
static Double_t GetBeautyMotherPtAndPDG(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Int_t &pdgGranma)
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 Bool_t IsTrackFromCharm(AliAODTrack *tr, TClonesArray *arrayMC)
static Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE)
Functions to check the decay tree.
TRandom * gRandom
static 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 Int_t CheckLbDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
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)
Double_t massD
static Double_t ResolK1(Double_t x)
static Int_t GetNumberOfTrackletsInEtaRange(AliAODEvent *ev, Double_t mineta, Double_t maxeta)
static Int_t CheckBsDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab, Bool_t ITS2UpgradeProd=kFALSE)
int Int_t
Definition: External.C:63
Double_t GetFullEvResol() const
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
static Int_t PreSelectITSUpgrade(TClonesArray *arrayMC, AliAODMCHeader *header, TObjArray aodTracks, Int_t nDaug, Int_t pdgabs, const Int_t *pdgDg)
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.
static Bool_t IsTrackFromHadronDecay(Int_t pdgMoth, AliAODTrack *tr, TClonesArray *arrayMC)
Definition: External.C:212
AliAODTrack * GetBachelor() const
static Double_t GetSubEvResolLowLim(const TH1F *hSubEvCorr)
const Double_t ptmin
UShort_t GetProngID(Int_t ip) const
static Bool_t HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header, TClonesArray *arrayMC)
static 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 Int_t CheckB0toDminuspiDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
static Double_t ComputeMaxd0MeasMinusExp(AliAODRecoDecayHF *cand, Double_t bfield)
Helper functions for D-meson analyses.
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)
static ROOT::Fit::FitResult DoInPlaneOutOfPlaneSimultaneusFit(AliHFInvMassFitter *&massfitterInPlane, AliHFInvMassFitter *&massfitterOutOfPlane, TH1F *hMassInPlane, TH1F *hMassOutOfPlane, Double_t MinMass, Double_t MaxMass, Double_t massD, vector< UInt_t > commonpars)
Double_t GetSubEvResol() const
static Bool_t IsTrackInjected(Int_t label, AliAODMCHeader *header, TClonesArray *arrayMC)
void SetInitialGaussianSigma(Double_t sigma)
void SetInitialGaussianMean(Double_t mean)
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 //.