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 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){
 244
 245   // Compute <pt> from 2D histogram M vs pt
 246
 247   //Make 2D histos in the desired pt range
 248   Int_t start=hMassD->FindBin(ptmin);
 249   Int_t end=hMassD->FindBin(ptmax)-1;
 250   const Int_t nx=end-start;
 251   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()));
 252   for(Int_t ix=start;ix<end;ix++){
 253     for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){
 254       hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy));
 255       hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy));
 256     }
 257   }
 258
 259   Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
 260   Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
 261   Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig);
 262   Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig);
 263   Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig);
 264   Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig);
 265   //  printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
 266
 267   Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
 268   Int_t minBinBkgLow=TMath::Max(hMassD->GetYaxis()->FindBin(minMass),2);
 269   Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
 270   Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
 271   Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
 272   Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
 273   Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
 274   Int_t maxBinBkgHi=TMath::Min(hMassD->GetYaxis()->FindBin(maxMass),hMassD->GetNbinsY()-1);
 275   Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
 276   Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
 277   //  printf("BKG Fit Limits = %f %f  && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
 278
 279   Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin);
 280   Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
 281   Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
 282   //  printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
 283
 284   TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow);
 285   TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi);
 286   TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig);
 287
 288   hMptBkgLo->Rebin(rebin);
 289   hMptBkgHi->Rebin(rebin);
 290   hMptSigReg->Rebin(rebin);
 291
 292   hMptBkgLo->Sumw2();
 293   hMptBkgHi->Sumw2();
 294   TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
 295   hMptBkgLoScal->Scale(bkgSig/bkgLow);
 296   TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
 297   hMptBkgHiScal->Scale(bkgSig/bkgHi);
 298
 299   TH1F* hMptBkgAver=0x0;
 300   hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin");
 301   hMptBkgAver->Add(hMptBkgHiScal);
 302   hMptBkgAver->Scale(0.5);
 303   TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin");
 304   hMptSig->Add(hMptBkgAver,-1.);
 305
 306   averagePt = hMptSig->GetMean();
 307   errorPt = hMptSig->GetMeanError();
 308
 309   delete hMptBkgLo;
 310   delete hMptBkgHi;
 311   delete hMptSigReg;
 312   delete hMptBkgLoScal;
 313   delete hMptBkgHiScal;
 314   delete hMptBkgAver;
 315   delete hMassDpt;
 316   delete hMptSig;
 317
 318 }
 319 //____________________________________________________________________________
 320 Bool_t AliVertexingHFUtils::CheckT0TriggerFired(AliAODEvent* aodEv){
 321   // check if T0VTX trigger was fired, based on a workaround suggested by Alla
 322   const Double32_t *mean = aodEv->GetT0TOF();
 323   if(mean && mean[0]<9999.) return kTRUE;
 324   else return kFALSE;
 325 }
 326 //____________________________________________________________________________
 327 Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
 328   // true impact parameter calculation for Dzero
 329
 330   if(!partD || !arrayMC || !mcHeader) return 99999.;
 331   Int_t code=partD->GetPdgCode();
 332   if(TMath::Abs(code)!=421) return 99999.;
 333
 334   Double_t vtxTrue[3];
 335   mcHeader->GetVertex(vtxTrue);
 336   Double_t origD[3];
 337   partD->XvYvZv(origD);
 338   Short_t charge=partD->Charge();
 339   Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
 340   for(Int_t iDau=0; iDau<2; iDau++){
 341     pXdauTrue[iDau]=0.;
 342     pYdauTrue[iDau]=0.;
 343     pZdauTrue[iDau]=0.;
 344   }
 345
 346   Int_t nDau=partD->GetNDaughters();
 347   Int_t labelFirstDau = partD->GetDaughter(0);
 348   if(nDau==2){
 349     for(Int_t iDau=0; iDau<2; iDau++){
 350       Int_t ind = labelFirstDau+iDau;
 351       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
 352       if(!part){
 353         printf("Daughter particle not found in MC array");
 354         return 99999.;
 355       }
 356       pXdauTrue[iDau]=part->Px();
 357       pYdauTrue[iDau]=part->Py();
 358       pZdauTrue[iDau]=part->Pz();
 359     }
 360   }else{
 361     return 99999.;
 362   }
 363
 364   Double_t d0dummy[2]={0.,0.};
 365   AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
 366   return aodDvsMC.ImpParXY();
 367
 368 }
 369 //____________________________________________________________________________
 370 Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
 371   // true impact parameter calculation for Dplus
 372
 373   if(!partD || !arrayMC || !mcHeader) return 99999.;
 374   Int_t code=partD->GetPdgCode();
 375   if(TMath::Abs(code)!=411) return 99999.;
 376
 377   Double_t vtxTrue[3];
 378   mcHeader->GetVertex(vtxTrue);
 379   Double_t origD[3];
 380   partD->XvYvZv(origD);
 381   Short_t charge=partD->Charge();
 382   Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
 383   for(Int_t iDau=0; iDau<3; iDau++){
 384     pXdauTrue[iDau]=0.;
 385     pYdauTrue[iDau]=0.;
 386     pZdauTrue[iDau]=0.;
 387   }
 388
 389   Int_t nDau=partD->GetNDaughters();
 390   Int_t labelFirstDau = partD->GetDaughter(0);
 391   if(nDau==3){
 392     for(Int_t iDau=0; iDau<3; iDau++){
 393       Int_t ind = labelFirstDau+iDau;
 394       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
 395       if(!part){
 396         printf("Daughter particle not found in MC array");
 397         return 99999.;
 398       }
 399       pXdauTrue[iDau]=part->Px();
 400       pYdauTrue[iDau]=part->Py();
 401       pZdauTrue[iDau]=part->Pz();
 402     }
 403   }else if(nDau==2){
 404     Int_t theDau=0;
 405     for(Int_t iDau=0; iDau<2; iDau++){
 406       Int_t ind = labelFirstDau+iDau;
 407       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
 408       if(!part){
 409         printf("Daughter particle not found in MC array");
 410         return 99999.;
 411       }
 412       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
 413       if(pdgCode==211 || pdgCode==321){
 414         pXdauTrue[theDau]=part->Px();
 415         pYdauTrue[theDau]=part->Py();
 416         pZdauTrue[theDau]=part->Pz();
 417         ++theDau;
 418       }else{
 419         Int_t nDauRes=part->GetNDaughters();
 420         if(nDauRes==2){
 421           Int_t labelFirstDauRes = part->GetDaughter(0);
 422           for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
 423             Int_t indDR = labelFirstDauRes+iDauRes;
 424             AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
 425             if(!partDR){
 426               printf("Daughter particle not found in MC array");
 427               return 99999.;
 428             }
 429
 430             Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
 431             if(pdgCodeDR==211 || pdgCodeDR==321){
 432               pXdauTrue[theDau]=partDR->Px();
 433               pYdauTrue[theDau]=partDR->Py();
 434               pZdauTrue[theDau]=partDR->Pz();
 435               ++theDau;
 436             }
 437           }
 438         }
 439       }
 440     }
 441     if(theDau!=3){
 442       printf("Wrong number of decay prongs");
 443       return 99999.;
 444     }
 445   }
 446
 447   Double_t d0dummy[3]={0.,0.,0.};
 448   AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
 449   return aodDvsMC.ImpParXY();
 450
 451 }
 452
 453
 454
 455 //____________________________________________________________________________
 456 Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) {
 457   //
 458   // Correct the number of accepted tracklets based on a TProfile Hist
 459   //
 460   //
 461
 462   if(TMath::Abs(vtxZ)>10.0){
 463     //    printf("ERROR: Z vertex out of range for correction of multiplicity\n");
 464     return uncorrectedNacc;
 465   }
 466
 467   if(!estimatorAvg){
 468     printf("ERROR: Missing TProfile for correction of multiplicity\n");
 469     return uncorrectedNacc;
 470   }
 471
 472   Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
 473
 474   Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
 475
 476   Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
 477
 478   return correctedNacc;
 479 }
 480 //______________________________________________________________________
 481 TString AliVertexingHFUtils::GetGenerator(Int_t label, AliAODMCHeader* header){
 482   // get the name of the generator that produced a given particle
 483
 484   Int_t nsumpart=0;
 485   TList *lh=header->GetCocktailHeaders();
 486   Int_t nh=lh->GetEntries();
 487   for(Int_t i=0;i<nh;i++){
 488     AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
 489     TString genname=gh->GetName();
 490     Int_t npart=gh->NProduced();
 491     if(label>=nsumpart && label<(nsumpart+npart)) return genname;
 492     nsumpart+=npart;
 493   }
 494   TString empty="";
 495   return empty;
 496 }
 497 //_____________________________________________________________________
 498 void AliVertexingHFUtils::GetTrackPrimaryGenerator(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
 499
 500   // method to check if a track comes from a given generator
 501
 502   Int_t lab=TMath::Abs(track->GetLabel());
 503   nameGen=GetGenerator(lab,header);
 504
 505   //  Int_t countControl=0;
 506
 507   while(nameGen.IsWhitespace()){
 508     AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
 509     if(!mcpart){
 510       printf("AliVertexingHFUtils::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
 511       break;
 512     }
 513     Int_t mother = mcpart->GetMother();
 514     if(mother<0){
 515       printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
 516       break;
 517     }
 518     lab=mother;
 519     nameGen=GetGenerator(mother,header);
 520     // countControl++;
 521     // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
 522     //   printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
 523     //   break;
 524     // }
 525   }
 526
 527   return;
 528 }
 529 //----------------------------------------------------------------------
 530 Bool_t AliVertexingHFUtils::IsTrackInjected(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC){
 531   // method to check if a track comes from the signal event or from the underlying Hijing event
 532   TString nameGen;
 533   GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
 534
 535   if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
 536
 537   return kTRUE;
 538 }
 539 //____________________________________________________________________________
 540 Bool_t AliVertexingHFUtils::IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
 541   // method to check if a D meson candidate comes from the signal event or from the underlying Hijing event
 542
 543   Int_t nprongs=cand->GetNProngs();
 544   for(Int_t i=0;i<nprongs;i++){
 545     AliAODTrack *daugh=(AliAODTrack*)cand->GetDaughter(i);
 546     if(IsTrackInjected(daugh,header,arrayMC)) return kTRUE;
 547   }
 548   return kFALSE;
 549 }
 550 //____________________________________________________________________________
 551 Bool_t AliVertexingHFUtils::HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
 552   // method to check if a cascade candidate comes from the signal event or from the underlying Hijing event
 553
 554   AliAODTrack* bach = cand->GetBachelor();
 555   if(IsTrackInjected(bach, header, arrayMC)) {
 556     AliDebug(2, "Bachelor is injected, the whole candidate is then injected");
 557     return kTRUE;
 558   }
 559   AliAODv0* v0 = cand->Getv0();
 560   Int_t nprongs = v0->GetNProngs();
 561   for(Int_t i = 0; i < nprongs; i++){
 562     AliAODTrack *daugh = (AliAODTrack*)v0->GetDaughter(i);
 563     if(IsTrackInjected(daugh,header,arrayMC)) {
 564       AliDebug(2, Form("V0 daughter number %d is injected, the whole candidate is then injected", i));
 565       return kTRUE;
 566     }
 567   }
 568   return kFALSE;
 569 }
 570 //____________________________________________________________________________
 571 Int_t AliVertexingHFUtils::CheckOrigin(AliStack* stack, TParticle *mcPart, Bool_t searchUpToQuark){
 572   // checking whether the mother of the particles come from a charm or a bottom quark
 573
 574   Int_t pdgGranma = 0;
 575   Int_t mother = 0;
 576   mother = mcPart->GetFirstMother();
 577   Int_t istep = 0;
 578   Int_t abspdgGranma =0;
 579   Bool_t isFromB=kFALSE;
 580   Bool_t isQuarkFound=kFALSE;
 581   while (mother >0 ){
 582     istep++;
 583     TParticle* mcGranma = stack->Particle(mother);
 584     if (mcGranma){
 585       pdgGranma = mcGranma->GetPdgCode();
 586       abspdgGranma = TMath::Abs(pdgGranma);
 587       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
 588         isFromB=kTRUE;
 589       }
 590       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
 591       mother = mcGranma->GetFirstMother();
 592     }else{
 593       printf("CheckOrigin: Failed casting the mother particle!");
 594       break;
 595     }
 596   }
 597   if(searchUpToQuark && !isQuarkFound) return 0;
 598   if(isFromB) return 5;
 599   else return 4;
 600
 601 }
 602 //____________________________________________________________________________
 603 Int_t AliVertexingHFUtils::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark){
 604   // checking whether the mother of the particles come from a charm or a bottom quark
 605
 606   Int_t pdgGranma = 0;
 607   Int_t mother = 0;
 608   mother = mcPart->GetMother();
 609   Int_t istep = 0;
 610   Int_t abspdgGranma =0;
 611   Bool_t isFromB=kFALSE;
 612   Bool_t isQuarkFound=kFALSE;
 613   while (mother >0 ){
 614     istep++;
 615     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
 616     if (mcGranma){
 617       pdgGranma = mcGranma->GetPdgCode();
 618       abspdgGranma = TMath::Abs(pdgGranma);
 619       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
 620         isFromB=kTRUE;
 621       }
 622       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
 623       mother = mcGranma->GetMother();
 624     }else{
 625       printf("AliVertexingHFUtils::CheckOrigin: Failed casting the mother particle!");
 626       break;
 627     }
 628   }
 629   if(searchUpToQuark && !isQuarkFound) return 0;
 630   if(isFromB) return 5;
 631   else return 4;
 632
 633 }
 634
 635 //____________________________________________________________________________
 636 Int_t AliVertexingHFUtils::CheckD0Decay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
 637   // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases
 638
 639   if(label<0) return -1;
 640   TParticle* mcPart = stack->Particle(label);
 641   Int_t pdgD=mcPart->GetPdgCode();
 642   if(TMath::Abs(pdgD)!=421) return -1;
 643
 644   Int_t nDau=mcPart->GetNDaughters();
 645
 646   if(nDau==2){
 647     Int_t daughter0 = mcPart->GetDaughter(0);
 648     Int_t daughter1 = mcPart->GetDaughter(1);
 649     TParticle* mcPartDaughter0 = stack->Particle(daughter0);
 650     TParticle* mcPartDaughter1 = stack->Particle(daughter1);
 651     if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
 652     arrayDauLab[0]=daughter0;
 653     arrayDauLab[1]=daughter1;
 654     Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
 655     Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
 656     if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
 657        !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
 658       return -1;
 659     }
 660     if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
 661     if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
 662     if((pdgCode0*pdgCode1)>0) return -1;
 663     Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
 664     Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
 665     Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
 666     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
 667     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
 668     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
 669     return 1;
 670   }
 671
 672   if(nDau==3 || nDau==4){
 673     Int_t nKaons=0;
 674     Int_t nPions=0;
 675     Double_t sumPxDau=0.;
 676     Double_t sumPyDau=0.;
 677     Double_t sumPzDau=0.;
 678     Int_t nFoundKpi=0;
 679     Int_t labelFirstDau = mcPart->GetDaughter(0);
 680     for(Int_t iDau=0; iDau<nDau; iDau++){
 681       Int_t indDau = labelFirstDau+iDau;
 682       if(indDau<0) return -1;
 683       TParticle* dau=stack->Particle(indDau);
 684       if(!dau) return -1;
 685       Int_t pdgdau=dau->GetPdgCode();
 686       if(TMath::Abs(pdgdau)==321){
 687         if(pdgD>0 && pdgdau>0) return -1;
 688         if(pdgD<0 && pdgdau<0) return -1;
 689         nKaons++;
 690         sumPxDau+=dau->Px();
 691         sumPyDau+=dau->Py();
 692         sumPzDau+=dau->Pz();
 693         arrayDauLab[nFoundKpi++]=indDau;
 694         if(nFoundKpi>4) return -1;
 695       }else if(TMath::Abs(pdgdau)==211){
 696         nPions++;
 697         sumPxDau+=dau->Px();
 698         sumPyDau+=dau->Py();
 699         sumPzDau+=dau->Pz();
 700         arrayDauLab[nFoundKpi++]=indDau;
 701         if(nFoundKpi>4) return -1;
 702       }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
 703         Int_t nResDau=dau->GetNDaughters();
 704         if(nResDau!=2) return -1;
 705         Int_t indFirstResDau=dau->GetDaughter(0);
 706         for(Int_t resDau=0; resDau<2; resDau++){
 707           Int_t indResDau=indFirstResDau+resDau;
 708           if(indResDau<0) return -1;
 709           TParticle* resdau=stack->Particle(indResDau);
 710           if(!resdau) return -1;
 711           Int_t pdgresdau=resdau->GetPdgCode();
 712           if(TMath::Abs(pdgresdau)==321){
 713             if(pdgD>0 && pdgresdau>0) return -1;
 714             if(pdgD<0 && pdgresdau<0) return -1;
 715             nKaons++;
 716             sumPxDau+=resdau->Px();
 717             sumPyDau+=resdau->Py();
 718             sumPzDau+=resdau->Pz();
 719             arrayDauLab[nFoundKpi++]=indResDau;
 720             if(nFoundKpi>4) return -1;
 721           }
 722           if(TMath::Abs(pdgresdau)==211){
 723             nPions++;
 724               sumPxDau+=resdau->Px();
 725               sumPyDau+=resdau->Py();
 726               sumPzDau+=resdau->Pz();
 727               arrayDauLab[nFoundKpi++]=indResDau;
 728               if(nFoundKpi>4) return -1;
 729           }
 730         }
 731       }else{
 732         return -1;
 733       }
 734     }
 735     if(nPions!=3) return -1;
 736     if(nKaons!=1) return -1;
 737     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
 738     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
 739     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
 740     return 2;
 741   }
 742   return -1;
 743 }
 744
 745 //____________________________________________________________________________
 746 Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
 747  // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases
 748
 749   Int_t pdgD=mcPart->GetPdgCode();
 750   if(TMath::Abs(pdgD)!=421) return -1;
 751
 752   Int_t nDau=mcPart->GetNDaughters();
 753
 754   if(nDau==2){
 755     Int_t daughter0 = mcPart->GetDaughter(0);
 756     Int_t daughter1 = mcPart->GetDaughter(1);
 757     AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter0));
 758     AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter1));
 759     if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
 760     arrayDauLab[0]=daughter0;
 761     arrayDauLab[1]=daughter1;
 762     Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
 763     Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
 764     if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
 765        !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
 766       return -1;
 767     }
 768     if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
 769     if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
 770     if((pdgCode0*pdgCode1)>0) return -1;
 771     Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
 772     Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
 773     Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
 774     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
 775     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
 776     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
 777     return 1;
 778   }
 779
 780   if(nDau==3 || nDau==4){
 781     Int_t nKaons=0;
 782     Int_t nPions=0;
 783     Double_t sumPxDau=0.;
 784     Double_t sumPyDau=0.;
 785     Double_t sumPzDau=0.;
 786     Int_t nFoundKpi=0;
 787     Int_t labelFirstDau = mcPart->GetDaughter(0);
 788     for(Int_t iDau=0; iDau<nDau; iDau++){
 789       Int_t indDau = labelFirstDau+iDau;
 790       if(indDau<0) return -1;
 791       AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
 792       if(!dau) return -1;
 793       Int_t pdgdau=dau->GetPdgCode();
 794       if(TMath::Abs(pdgdau)==321){
 795         if(pdgD>0 && pdgdau>0) return -1;
 796         if(pdgD<0 && pdgdau<0) return -1;
 797         nKaons++;
 798         sumPxDau+=dau->Px();
 799         sumPyDau+=dau->Py();
 800         sumPzDau+=dau->Pz();
 801         arrayDauLab[nFoundKpi++]=indDau;
 802         if(nFoundKpi>4) return -1;
 803       }else if(TMath::Abs(pdgdau)==211){
 804         nPions++;
 805         sumPxDau+=dau->Px();
 806         sumPyDau+=dau->Py();
 807         sumPzDau+=dau->Pz();
 808         arrayDauLab[nFoundKpi++]=indDau;
 809         if(nFoundKpi>4) return -1;
 810       }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
 811         Int_t nResDau=dau->GetNDaughters();
 812         if(nResDau!=2) return -1;
 813         Int_t indFirstResDau=dau->GetDaughter(0);
 814         for(Int_t resDau=0; resDau<2; resDau++){
 815           Int_t indResDau=indFirstResDau+resDau;
 816           if(indResDau<0) return -1;
 817           AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
 818           if(!resdau) return -1;
 819           Int_t pdgresdau=resdau->GetPdgCode();
 820           if(TMath::Abs(pdgresdau)==321){
 821             if(pdgD>0 && pdgresdau>0) return -1;
 822             if(pdgD<0 && pdgresdau<0) return -1;
 823             nKaons++;
 824             sumPxDau+=resdau->Px();
 825             sumPyDau+=resdau->Py();
 826             sumPzDau+=resdau->Pz();
 827             arrayDauLab[nFoundKpi++]=indResDau;
 828             if(nFoundKpi>4) return -1;
 829           }
 830           if(TMath::Abs(pdgresdau)==211){
 831             nPions++;
 832             sumPxDau+=resdau->Px();
 833             sumPyDau+=resdau->Py();
 834             sumPzDau+=resdau->Pz();
 835             arrayDauLab[nFoundKpi++]=indResDau;
 836             if(nFoundKpi>4) return -1;
 837           }
 838         }
 839       }else{
 840         return -1;
 841       }
 842     }
 843     if(nPions!=3) return -1;
 844     if(nKaons!=1) return -1;
 845     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
 846     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
 847     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
 848     return 2;
 849   }
 850   return -1;
 851 }
 852 //____________________________________________________________________________
 853 Int_t AliVertexingHFUtils::CheckDplusDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
 854   // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases
 855
 856   if(label<0) return -1;
 857   TParticle* mcPart = stack->Particle(label);
 858   Int_t pdgD=mcPart->GetPdgCode();
 859   if(TMath::Abs(pdgD)!=411) return -1;
 860
 861   Int_t nDau=mcPart->GetNDaughters();
 862   Int_t labelFirstDau = mcPart->GetDaughter(0);
 863   Int_t nKaons=0;
 864   Int_t nPions=0;
 865   Double_t sumPxDau=0.;
 866   Double_t sumPyDau=0.;
 867   Double_t sumPzDau=0.;
 868   Int_t nFoundKpi=0;
 869
 870   if(nDau==3 || nDau==2){
 871     for(Int_t iDau=0; iDau<nDau; iDau++){
 872       Int_t indDau = labelFirstDau+iDau;
 873       if(indDau<0) return -1;
 874       TParticle* dau=stack->Particle(indDau);
 875       if(!dau) return -1;
 876       Int_t pdgdau=dau->GetPdgCode();
 877       if(TMath::Abs(pdgdau)==321){
 878         if(pdgD*pdgdau>0) return -1;
 879         nKaons++;
 880         sumPxDau+=dau->Px();
 881         sumPyDau+=dau->Py();
 882         sumPzDau+=dau->Pz();
 883         arrayDauLab[nFoundKpi++]=indDau;
 884         if(nFoundKpi>3) return -1;
 885       }else if(TMath::Abs(pdgdau)==211){
 886         if(pdgD*pdgdau<0) return -1;
 887         nPions++;
 888         sumPxDau+=dau->Px();
 889         sumPyDau+=dau->Py();
 890         sumPzDau+=dau->Pz();
 891         arrayDauLab[nFoundKpi++]=indDau;
 892         if(nFoundKpi>3) return -1;
 893       }else if(TMath::Abs(pdgdau)==313){
 894         Int_t nResDau=dau->GetNDaughters();
 895         if(nResDau!=2) return -1;
 896         Int_t indFirstResDau=dau->GetDaughter(0);
 897         for(Int_t resDau=0; resDau<2; resDau++){
 898           Int_t indResDau=indFirstResDau+resDau;
 899           if(indResDau<0) return -1;
 900           TParticle* resdau=stack->Particle(indResDau);
 901           if(!resdau) return -1;
 902           Int_t pdgresdau=resdau->GetPdgCode();
 903           if(TMath::Abs(pdgresdau)==321){
 904             if(pdgD*pdgresdau>0) return -1;
 905             sumPxDau+=resdau->Px();
 906             sumPyDau+=resdau->Py();
 907             sumPzDau+=resdau->Pz();
 908             nKaons++;
 909             arrayDauLab[nFoundKpi++]=indResDau;
 910             if(nFoundKpi>3) return -1;
 911           }
 912           if(TMath::Abs(pdgresdau)==211){
 913             if(pdgD*pdgresdau<0) return -1;
 914             sumPxDau+=resdau->Px();
 915             sumPyDau+=resdau->Py();
 916             sumPzDau+=resdau->Pz();
 917             nPions++;
 918             arrayDauLab[nFoundKpi++]=indResDau;
 919             if(nFoundKpi>3) return -1;
 920           }
 921           }
 922       }else{
 923         return -1;
 924       }
 925     }
 926     if(nPions!=2) return -1;
 927     if(nKaons!=1) return -1;
 928     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
 929     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
 930     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
 931     if(nDau==3) return 1;
 932     else if(nDau==2) return 2;
 933   }
 934
 935   return -1;
 936
 937 }
 938
 939 //____________________________________________________________________________
 940 Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
 941   // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases
 942
 943   Int_t pdgD=mcPart->GetPdgCode();
 944   if(TMath::Abs(pdgD)!=411) return -1;
 945
 946   Int_t nDau=mcPart->GetNDaughters();
 947   Int_t labelFirstDau = mcPart->GetDaughter(0);
 948   Int_t nKaons=0;
 949   Int_t nPions=0;
 950   Double_t sumPxDau=0.;
 951   Double_t sumPyDau=0.;
 952   Double_t sumPzDau=0.;
 953   Int_t nFoundKpi=0;
 954
 955   if(nDau==3 || nDau==2){
 956     for(Int_t iDau=0; iDau<nDau; iDau++){
 957       Int_t indDau = labelFirstDau+iDau;
 958       if(indDau<0) return -1;
 959       AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
 960       if(!dau) return -1;
 961       Int_t pdgdau=dau->GetPdgCode();
 962       if(TMath::Abs(pdgdau)==321){
 963         if(pdgD*pdgdau>0) return -1;
 964         nKaons++;
 965         sumPxDau+=dau->Px();
 966         sumPyDau+=dau->Py();
 967         sumPzDau+=dau->Pz();
 968         arrayDauLab[nFoundKpi++]=indDau;
 969         if(nFoundKpi>3) return -1;
 970       }else if(TMath::Abs(pdgdau)==211){
 971         if(pdgD*pdgdau<0) return -1;
 972         nPions++;
 973         sumPxDau+=dau->Px();
 974         sumPyDau+=dau->Py();
 975         sumPzDau+=dau->Pz();
 976         arrayDauLab[nFoundKpi++]=indDau;
 977         if(nFoundKpi>3) return -1;
 978       }else if(TMath::Abs(pdgdau)==313){
 979         Int_t nResDau=dau->GetNDaughters();
 980         if(nResDau!=2) return -1;
 981         Int_t indFirstResDau=dau->GetDaughter(0);
 982         for(Int_t resDau=0; resDau<2; resDau++){
 983           Int_t indResDau=indFirstResDau+resDau;
 984           if(indResDau<0) return -1;
 985           AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
 986           if(!resdau) return -1;
 987           Int_t pdgresdau=resdau->GetPdgCode();
 988           if(TMath::Abs(pdgresdau)==321){
 989             if(pdgD*pdgresdau>0) return -1;
 990             sumPxDau+=resdau->Px();
 991             sumPyDau+=resdau->Py();
 992             sumPzDau+=resdau->Pz();
 993             nKaons++;
 994             arrayDauLab[nFoundKpi++]=indResDau;
 995             if(nFoundKpi>3) return -1;
 996           }
 997           if(TMath::Abs(pdgresdau)==211){
 998             if(pdgD*pdgresdau<0) return -1;
 999             sumPxDau+=resdau->Px();
1000             sumPyDau+=resdau->Py();
1001             sumPzDau+=resdau->Pz();
1002             nPions++;
1003             arrayDauLab[nFoundKpi++]=indResDau;
1004             if(nFoundKpi>3) return -1;
1005           }
1006         }
1007        }else{
1008         return -1;
1009       }
1010     }
1011     if(nPions!=2) return -1;
1012     if(nKaons!=1) return -1;
1013     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1014     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1015     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1016     if(nDau==3) return 1;
1017     else if(nDau==2) return 2;
1018   }
1019
1020   return -1;
1021
1022 }
1023 //____________________________________________________________________________
1024 Int_t AliVertexingHFUtils::CheckDsDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1025   // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case
1026
1027   if(label<0) return -1;
1028   TParticle* mcPart = stack->Particle(label);
1029   Int_t pdgD=mcPart->GetPdgCode();
1030   if(TMath::Abs(pdgD)!=431) return -1;
1031
1032   Int_t nDau=mcPart->GetNDaughters();
1033   Int_t labelFirstDau = mcPart->GetDaughter(0);
1034   Int_t nKaons=0;
1035   Int_t nPions=0;
1036   Double_t sumPxDau=0.;
1037   Double_t sumPyDau=0.;
1038   Double_t sumPzDau=0.;
1039   Int_t nFoundKpi=0;
1040   Bool_t isPhi=kFALSE;
1041   Bool_t isk0st=kFALSE;
1042
1043   if(nDau==3 || nDau==2){
1044     for(Int_t iDau=0; iDau<nDau; iDau++){
1045       Int_t indDau = labelFirstDau+iDau;
1046       if(indDau<0) return -1;
1047       TParticle* dau=stack->Particle(indDau);
1048       if(!dau) return -1;
1049       Int_t pdgdau=dau->GetPdgCode();
1050       if(TMath::Abs(pdgdau)==321){
1051         nKaons++;
1052         sumPxDau+=dau->Px();
1053         sumPyDau+=dau->Py();
1054         sumPzDau+=dau->Pz();
1055         arrayDauLab[nFoundKpi++]=indDau;
1056         if(nFoundKpi>3) return -1;
1057       }else if(TMath::Abs(pdgdau)==211){
1058         nPions++;
1059         sumPxDau+=dau->Px();
1060         sumPyDau+=dau->Py();
1061         sumPzDau+=dau->Pz();
1062         arrayDauLab[nFoundKpi++]=indDau;
1063         if(nFoundKpi>3) return -1;
1064       }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1065         if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1066         if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1067         Int_t nResDau=dau->GetNDaughters();
1068         if(nResDau!=2) return -1;
1069         Int_t indFirstResDau=dau->GetDaughter(0);
1070         for(Int_t resDau=0; resDau<2; resDau++){
1071           Int_t indResDau=indFirstResDau+resDau;
1072           if(indResDau<0) return -1;
1073           TParticle* resdau=stack->Particle(indResDau);
1074           if(!resdau) return -1;
1075           Int_t pdgresdau=resdau->GetPdgCode();
1076           if(TMath::Abs(pdgresdau)==321){
1077             sumPxDau+=resdau->Px();
1078             sumPyDau+=resdau->Py();
1079             sumPzDau+=resdau->Pz();
1080             nKaons++;
1081             arrayDauLab[nFoundKpi++]=indResDau;
1082             if(nFoundKpi>3) return -1;
1083           }
1084           if(TMath::Abs(pdgresdau)==211){
1085             sumPxDau+=resdau->Px();
1086             sumPyDau+=resdau->Py();
1087             sumPzDau+=resdau->Pz();
1088             nPions++;
1089             arrayDauLab[nFoundKpi++]=indResDau;
1090             if(nFoundKpi>3) return -1;
1091           }
1092           }
1093       }else{
1094         return -1;
1095       }
1096     }
1097     if(nPions!=1) return -1;
1098     if(nKaons!=2) return -1;
1099     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1100     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1101     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1102     if(nDau==3) return 3;
1103     else if(nDau==2){
1104       if(isk0st) return 2;
1105       if(isPhi) return 1;
1106     }
1107   }
1108
1109   return -1;
1110
1111 }
1112 //____________________________________________________________________________
1113 Int_t AliVertexingHFUtils::CheckDsDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1114   // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case
1115
1116   Int_t pdgD=mcPart->GetPdgCode();
1117   if(TMath::Abs(pdgD)!=431) return -1;
1118
1119   Int_t nDau=mcPart->GetNDaughters();
1120   Int_t labelFirstDau = mcPart->GetDaughter(0);
1121   Int_t nKaons=0;
1122   Int_t nPions=0;
1123   Double_t sumPxDau=0.;
1124   Double_t sumPyDau=0.;
1125   Double_t sumPzDau=0.;
1126   Int_t nFoundKpi=0;
1127   Bool_t isPhi=kFALSE;
1128   Bool_t isk0st=kFALSE;
1129
1130   if(nDau==3 || nDau==2){
1131     for(Int_t iDau=0; iDau<nDau; iDau++){
1132       Int_t indDau = labelFirstDau+iDau;
1133       if(indDau<0) return -1;
1134       AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1135       if(!dau) return -1;
1136       Int_t pdgdau=dau->GetPdgCode();
1137       if(TMath::Abs(pdgdau)==321){
1138         nKaons++;
1139         sumPxDau+=dau->Px();
1140         sumPyDau+=dau->Py();
1141         sumPzDau+=dau->Pz();
1142         arrayDauLab[nFoundKpi++]=indDau;
1143         if(nFoundKpi>3) return -1;
1144       }else if(TMath::Abs(pdgdau)==211){
1145         nPions++;
1146         sumPxDau+=dau->Px();
1147         sumPyDau+=dau->Py();
1148         sumPzDau+=dau->Pz();
1149         arrayDauLab[nFoundKpi++]=indDau;
1150         if(nFoundKpi>3) return -1;
1151       }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1152         if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1153         if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1154         Int_t nResDau=dau->GetNDaughters();
1155         if(nResDau!=2) return -1;
1156         Int_t indFirstResDau=dau->GetDaughter(0);
1157         for(Int_t resDau=0; resDau<2; resDau++){
1158           Int_t indResDau=indFirstResDau+resDau;
1159           if(indResDau<0) return -1;
1160           AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1161           if(!resdau) return -1;
1162           Int_t pdgresdau=resdau->GetPdgCode();
1163           if(TMath::Abs(pdgresdau)==321){
1164             sumPxDau+=resdau->Px();
1165             sumPyDau+=resdau->Py();
1166             sumPzDau+=resdau->Pz();
1167             nKaons++;
1168             arrayDauLab[nFoundKpi++]=indResDau;
1169             if(nFoundKpi>3) return -1;
1170           }
1171           if(TMath::Abs(pdgresdau)==211){
1172             sumPxDau+=resdau->Px();
1173             sumPyDau+=resdau->Py();
1174             sumPzDau+=resdau->Pz();
1175             nPions++;
1176             arrayDauLab[nFoundKpi++]=indResDau;
1177             if(nFoundKpi>3) return -1;
1178           }
1179           }
1180       }else{
1181         return -1;
1182       }
1183     }
1184     if(nPions!=1) return -1;
1185     if(nKaons!=2) return -1;
1186     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1187     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1188     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1189     if(nDau==3) return 3;
1190     else if(nDau==2){
1191       if(isk0st) return 2;
1192       if(isPhi) return 1;
1193     }
1194   }
1195
1196   return -1;
1197
1198 }
1199 //____________________________________________________________________________
1200 Int_t AliVertexingHFUtils::CheckDstarDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1201   // Checks the Dstar decay channel. Returns 1 for D*->D0pi->Kpipi, -1 in other cases
1202
1203   if(label<0) return -1;
1204   TParticle* mcPart = stack->Particle(label);
1205   Int_t pdgD=mcPart->GetPdgCode();
1206   if(TMath::Abs(pdgD)!=413) return -1;
1207
1208   Int_t nDau=mcPart->GetNDaughters();
1209   if(nDau!=2) return -1;
1210
1211   Int_t labelFirstDau = mcPart->GetDaughter(0);
1212   Int_t nKaons=0;
1213   Int_t nPions=0;
1214   Double_t sumPxDau=0.;
1215   Double_t sumPyDau=0.;
1216   Double_t sumPzDau=0.;
1217   Int_t nFoundKpi=0;
1218
1219   for(Int_t iDau=0; iDau<nDau; iDau++){
1220     Int_t indDau = labelFirstDau+iDau;
1221     if(indDau<0) return -1;
1222     TParticle* dau=stack->Particle(indDau);
1223     if(!dau) return -1;
1224     Int_t pdgdau=dau->GetPdgCode();
1225     if(TMath::Abs(pdgdau)==421){
1226       Int_t nResDau=dau->GetNDaughters();
1227       if(nResDau!=2) return -1;
1228       Int_t indFirstResDau=dau->GetDaughter(0);
1229       for(Int_t resDau=0; resDau<2; resDau++){
1230         Int_t indResDau=indFirstResDau+resDau;
1231         if(indResDau<0) return -1;
1232         TParticle* resdau=stack->Particle(indResDau);
1233         if(!resdau) return -1;
1234         Int_t pdgresdau=resdau->GetPdgCode();
1235         if(TMath::Abs(pdgresdau)==321){
1236           if(pdgD*pdgresdau>0) return -1;
1237           sumPxDau+=resdau->Px();
1238           sumPyDau+=resdau->Py();
1239           sumPzDau+=resdau->Pz();
1240           nKaons++;
1241           arrayDauLab[nFoundKpi++]=indResDau;
1242           if(nFoundKpi>3) return -1;
1243         }
1244         if(TMath::Abs(pdgresdau)==211){
1245           if(pdgD*pdgresdau<0) return -1;
1246           sumPxDau+=resdau->Px();
1247           sumPyDau+=resdau->Py();
1248           sumPzDau+=resdau->Pz();
1249           nPions++;
1250           arrayDauLab[nFoundKpi++]=indResDau;
1251           if(nFoundKpi>3) return -1;
1252         }
1253       }
1254     }else if(TMath::Abs(pdgdau)==211){
1255       if(pdgD*pdgdau<0) return -1;
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     }
1263   }
1264
1265   if(nPions!=2) return -1;
1266   if(nKaons!=1) return -1;
1267   if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1268   if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1269   if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1270   return 1;
1271
1272 }
1273
1274 //____________________________________________________________________________
1275 Int_t AliVertexingHFUtils::CheckDstarDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1276   // Checks the Dstar decay channel. Returns 1 for D*->D0pi->Kpipi, -1 in other cases
1277
1278   Int_t pdgD=mcPart->GetPdgCode();
1279   if(TMath::Abs(pdgD)!=413) return -1;
1280
1281   Int_t nDau=mcPart->GetNDaughters();
1282   if(nDau!=2) return -1;
1283
1284   Int_t labelFirstDau = mcPart->GetDaughter(0);
1285   Int_t nKaons=0;
1286   Int_t nPions=0;
1287   Double_t sumPxDau=0.;
1288   Double_t sumPyDau=0.;
1289   Double_t sumPzDau=0.;
1290   Int_t nFoundKpi=0;
1291
1292   for(Int_t iDau=0; iDau<nDau; iDau++){
1293     Int_t indDau = labelFirstDau+iDau;
1294     if(indDau<0) return -1;
1295     AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1296     if(!dau) return -1;
1297     Int_t pdgdau=dau->GetPdgCode();
1298     if(TMath::Abs(pdgdau)==421){
1299       Int_t nResDau=dau->GetNDaughters();
1300       if(nResDau!=2) return -1;
1301       Int_t indFirstResDau=dau->GetDaughter(0);
1302       for(Int_t resDau=0; resDau<2; resDau++){
1303         Int_t indResDau=indFirstResDau+resDau;
1304         if(indResDau<0) return -1;
1305         AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1306         if(!resdau) return -1;
1307         Int_t pdgresdau=resdau->GetPdgCode();
1308         if(TMath::Abs(pdgresdau)==321){
1309           if(pdgD*pdgresdau>0) return -1;
1310           sumPxDau+=resdau->Px();
1311           sumPyDau+=resdau->Py();
1312           sumPzDau+=resdau->Pz();
1313           nKaons++;
1314           arrayDauLab[nFoundKpi++]=indResDau;
1315           if(nFoundKpi>3) return -1;
1316         }
1317         if(TMath::Abs(pdgresdau)==211){
1318           if(pdgD*pdgresdau<0) return -1;
1319           sumPxDau+=resdau->Px();
1320           sumPyDau+=resdau->Py();
1321           sumPzDau+=resdau->Pz();
1322           nPions++;
1323           arrayDauLab[nFoundKpi++]=indResDau;
1324           if(nFoundKpi>3) return -1;
1325         }
1326       }
1327     }else if(TMath::Abs(pdgdau)==211){
1328       if(pdgD*pdgdau<0) return -1;
1329       nPions++;
1330       sumPxDau+=dau->Px();
1331       sumPyDau+=dau->Py();
1332       sumPzDau+=dau->Pz();
1333       arrayDauLab[nFoundKpi++]=indDau;
1334       if(nFoundKpi>3) return -1;
1335     }
1336   }
1337
1338   if(nPions!=2) return -1;
1339   if(nKaons!=1) return -1;
1340   if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1341   if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1342   if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1343   return 1;
1344
1345 }
1346