PWGHF/vertexingHF/AliVertexingHFUtils.cxx

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