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 "AliAODEvent.h"
  24 #include "AliAODMCHeader.h"
  25 #include "AliGenEventHeader.h"
  26 #include "AliAODMCParticle.h"
  27 #include "AliAODRecoDecayHF.h"
  28 #include "AliVertexingHFUtils.h"
  29
  30 /* $Id$ */
  31
  32 ///////////////////////////////////////////////////////////////////
  33 //                                                               //
  34 // Class with functions useful for different D2H analyses        //
  35 // - event plane resolution                                      //
  36 // - <pt> calculation with side band subtraction                 //
  37 // - tracklet multiplicity calculation                            //
  38 // Origin: F.Prino, Torino, prino@to.infn.it                     //
  39 //                                                               //
  40 ///////////////////////////////////////////////////////////////////
  41
  42 ClassImp(AliVertexingHFUtils)
  43
  44 //______________________________________________________________________
  45 AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
  46   fK(1),
  47   fSubRes(1.),
  48   fMinEtaForTracklets(-1.),
  49   fMaxEtaForTracklets(1.)
  50 {
  51   // Default contructor
  52 }
  53
  54
  55 //______________________________________________________________________
  56 AliVertexingHFUtils::AliVertexingHFUtils(Int_t k):
  57   TObject(),
  58   fK(k),
  59   fSubRes(1.),
  60   fMinEtaForTracklets(-1.),
  61   fMaxEtaForTracklets(1.)
  62 {
  63   // Standard constructor
  64 }
  65
  66
  67 //______________________________________________________________________
  68 void AliVertexingHFUtils::ComputeSignificance(Double_t signal, Double_t  errsignal, Double_t  background, Double_t  errbackground, Double_t &significance,Double_t &errsignificance){
  69   // calculate significance from S, B and errors
  70
  71
  72   Double_t errSigSq=errsignal*errsignal;
  73   Double_t errBkgSq=errbackground*errbackground;
  74   Double_t sigPlusBkg=signal+background;
  75   if(sigPlusBkg>0. && signal>0.){
  76     significance =  signal/TMath::Sqrt(signal+background);
  77     errsignificance = significance*TMath::Sqrt((errSigSq+errBkgSq)/(4.*sigPlusBkg*sigPlusBkg)+(background/sigPlusBkg)*errSigSq/signal/signal);
  78   }else{
  79     significance=0.;
  80     errsignificance=0.;
  81   }
  82   return;
  83
  84 }
  85 //______________________________________________________________________
  86 Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){
  87   // compute chi from polynomial approximation
  88   Double_t c[5];
  89   if(k==1){
  90     c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
  91   }
  92   else if(k==2){
  93     c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
  94   }
  95   else return -1;
  96   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;
  97 }
  98
  99 //______________________________________________________________________
 100 Double_t AliVertexingHFUtils:: ResolK1(Double_t x){
 101   return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
 102 }
 103
 104
 105 //______________________________________________________________________
 106 Double_t AliVertexingHFUtils::FindChi(Double_t res,  Int_t k){
 107   // compute chi variable (=v2*sqrt(N)) from external values
 108
 109   Double_t x1=0;
 110   Double_t x2=15;
 111   Double_t y1,y2;
 112   if(k==1){
 113     y1=ResolK1(x1)-res;
 114     y2=ResolK1(x2)-res;
 115   }
 116   else if(k==2){
 117     y1=Pol(x1,2)-res;
 118     y2=Pol(x2,2)-res;
 119   }
 120   else return -1;
 121
 122   if(y1*y2>0) return -1;
 123   if(y1==0) return y1;
 124   if(y2==0) return y2;
 125   Double_t xmed,ymed;
 126   Int_t jiter=0;
 127   while((x2-x1)>0.0001){
 128     xmed=0.5*(x1+x2);
 129     if(k==1){
 130       y1=ResolK1(x1)-res;
 131       y2=ResolK1(x2)-res;
 132       ymed=ResolK1(xmed)-res;
 133     }
 134     else if(k==2){
 135       y1=Pol(x1,2)-res;
 136       y2=Pol(x2,2)-res;
 137       ymed=Pol(xmed,2)-res;
 138     }
 139     else return -1;
 140     if((y1*ymed)<0) x2=xmed;
 141     if((y2*ymed)<0)x1=xmed;
 142     if(ymed==0) return xmed;
 143     jiter++;
 144   }
 145   return 0.5*(x1+x2);
 146 }
 147
 148 //______________________________________________________________________
 149 Double_t AliVertexingHFUtils::GetFullEvResol(Double_t resSub, Int_t k){
 150   // computes event plane resolution starting from sub event resolution
 151   Double_t chisub=FindChi(resSub,k);
 152   Double_t chifull=chisub*TMath::Sqrt(2);
 153   if(k==1) return ResolK1(chifull);
 154   else if(k==2) return Pol(chifull,2);
 155   else{
 156     printf("k should be <=2\n");
 157     return 1.;
 158   }
 159 }
 160
 161 //______________________________________________________________________
 162 Double_t AliVertexingHFUtils::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){
 163   // computes event plane resolution starting from sub event correlation histogram
 164   if(!hSubEvCorr) return 1.;
 165   Double_t resSub=GetSubEvResol(hSubEvCorr);
 166   return GetFullEvResol(resSub,k);
 167 }
 168 //______________________________________________________________________
 169 Double_t AliVertexingHFUtils::GetFullEvResolLowLim(const TH1F* hSubEvCorr, Int_t k){
 170   // computes low limit event plane resolution starting from sub event correlation histogram
 171   if(!hSubEvCorr) return 1.;
 172   Double_t resSub=GetSubEvResolLowLim(hSubEvCorr);
 173   printf("%f\n",resSub);
 174   return GetFullEvResol(resSub,k);
 175 }
 176 //______________________________________________________________________
 177 Double_t AliVertexingHFUtils::GetFullEvResolHighLim(const TH1F* hSubEvCorr, Int_t k){
 178   // computes high limit event plane resolution starting from sub event correlation histogram
 179   if(!hSubEvCorr) return 1.;
 180   Double_t resSub=GetSubEvResolHighLim(hSubEvCorr);
 181   printf("%f\n",resSub);
 182   return GetFullEvResol(resSub,k);
 183 }
 184 //______________________________________________________________________
 185 Int_t AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(AliAODEvent* ev, Double_t mineta, Double_t maxeta){
 186   // counts tracklets in given eta range
 187   AliAODTracklets* tracklets=ev->GetTracklets();
 188   Int_t nTr=tracklets->GetNumberOfTracklets();
 189   Int_t count=0;
 190   for(Int_t iTr=0; iTr<nTr; iTr++){
 191     Double_t theta=tracklets->GetTheta(iTr);
 192     Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
 193     if(eta>mineta && eta<maxeta) count++;
 194   }
 195   return count;
 196 }
 197 //______________________________________________________________________
 198 Int_t AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
 199   // counts generated particles in fgiven eta range
 200
 201   Int_t nChargedMC=0;
 202   for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
 203     AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
 204     Int_t charge = part->Charge();
 205     Double_t eta = part->Eta();
 206     if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++;
 207   }
 208   return nChargedMC;
 209 }
 210 //______________________________________________________________________
 211 Int_t AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
 212   // counts generated primary particles in given eta range
 213
 214   Int_t nChargedMC=0;
 215   for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
 216     AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
 217     Int_t charge = part->Charge();
 218     Double_t eta = part->Eta();
 219     if(charge!=0 && eta>mineta && eta<maxeta){
 220       if(part->IsPrimary())nChargedMC++;
 221     }
 222   }
 223   return nChargedMC;
 224 }
 225 //______________________________________________________________________
 226 Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
 227   // counts generated primary particles in given eta range
 228
 229   Int_t nChargedMC=0;
 230   for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
 231     AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
 232     Int_t charge = part->Charge();
 233     Double_t eta = part->Eta();
 234     if(charge!=0 && eta>mineta && eta<maxeta){
 235       if(part->IsPhysicalPrimary())nChargedMC++;
 236     }
 237   }
 238   return nChargedMC;
 239 }
 240 //______________________________________________________________________
 241 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){
 242
 243   // Compute <pt> from 2D histogram M vs pt
 244
 245   //Make 2D histos in the desired pt range
 246   Int_t start=hMassD->FindBin(ptmin);
 247   Int_t end=hMassD->FindBin(ptmax)-1;
 248   const Int_t nx=end-start;
 249   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()));
 250   for(Int_t ix=start;ix<end;ix++){
 251     for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){
 252       hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy));
 253       hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy));
 254     }
 255   }
 256
 257   Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
 258   Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
 259   Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig);
 260   Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig);
 261   Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig);
 262   Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig);
 263   //  printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
 264
 265   Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
 266   Int_t minBinBkgLow=TMath::Max(hMassD->GetYaxis()->FindBin(minMass),2);
 267   Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
 268   Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
 269   Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
 270   Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
 271   Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
 272   Int_t maxBinBkgHi=TMath::Min(hMassD->GetYaxis()->FindBin(maxMass),hMassD->GetNbinsY()-1);
 273   Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
 274   Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
 275   //  printf("BKG Fit Limits = %f %f  && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
 276
 277   Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin);
 278   Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
 279   Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
 280   //  printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
 281
 282   TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow);
 283   TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi);
 284   TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig);
 285
 286   hMptBkgLo->Rebin(rebin);
 287   hMptBkgHi->Rebin(rebin);
 288   hMptSigReg->Rebin(rebin);
 289
 290   hMptBkgLo->Sumw2();
 291   hMptBkgHi->Sumw2();
 292   TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
 293   hMptBkgLoScal->Scale(bkgSig/bkgLow);
 294   TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
 295   hMptBkgHiScal->Scale(bkgSig/bkgHi);
 296
 297   TH1F* hMptBkgAver=0x0;
 298   hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin");
 299   hMptBkgAver->Add(hMptBkgHiScal);
 300   hMptBkgAver->Scale(0.5);
 301   TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin");
 302   hMptSig->Add(hMptBkgAver,-1.);
 303
 304   averagePt = hMptSig->GetMean();
 305   errorPt = hMptSig->GetMeanError();
 306
 307   delete hMptBkgLo;
 308   delete hMptBkgHi;
 309   delete hMptSigReg;
 310   delete hMptBkgLoScal;
 311   delete hMptBkgHiScal;
 312   delete hMptBkgAver;
 313   delete hMassDpt;
 314   delete hMptSig;
 315
 316 }
 317 //____________________________________________________________________________
 318 Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
 319   // true impact parameter calculation for Dzero
 320
 321   if(!partD || !arrayMC || !mcHeader) return 99999.;
 322   Int_t code=partD->GetPdgCode();
 323   if(TMath::Abs(code)!=421) return 99999.;
 324
 325   Double_t vtxTrue[3];
 326   mcHeader->GetVertex(vtxTrue);
 327   Double_t origD[3];
 328   partD->XvYvZv(origD);
 329   Short_t charge=partD->Charge();
 330   Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
 331   for(Int_t iDau=0; iDau<2; iDau++){
 332     pXdauTrue[iDau]=0.;
 333     pYdauTrue[iDau]=0.;
 334     pZdauTrue[iDau]=0.;
 335   }
 336
 337   Int_t nDau=partD->GetNDaughters();
 338   Int_t labelFirstDau = partD->GetDaughter(0);
 339   if(nDau==2){
 340     for(Int_t iDau=0; iDau<2; iDau++){
 341       Int_t ind = labelFirstDau+iDau;
 342       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
 343       if(!part){
 344         printf("Daughter particle not found in MC array");
 345         return 99999.;
 346       }
 347       pXdauTrue[iDau]=part->Px();
 348       pYdauTrue[iDau]=part->Py();
 349       pZdauTrue[iDau]=part->Pz();
 350     }
 351   }else{
 352     return 99999.;
 353   }
 354
 355   Double_t d0dummy[2]={0.,0.};
 356   AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
 357   return aodDvsMC.ImpParXY();
 358
 359 }
 360 //____________________________________________________________________________
 361 Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
 362   // true impact parameter calculation for Dplus
 363
 364   if(!partD || !arrayMC || !mcHeader) return 99999.;
 365   Int_t code=partD->GetPdgCode();
 366   if(TMath::Abs(code)!=411) return 99999.;
 367
 368   Double_t vtxTrue[3];
 369   mcHeader->GetVertex(vtxTrue);
 370   Double_t origD[3];
 371   partD->XvYvZv(origD);
 372   Short_t charge=partD->Charge();
 373   Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
 374   for(Int_t iDau=0; iDau<3; iDau++){
 375     pXdauTrue[iDau]=0.;
 376     pYdauTrue[iDau]=0.;
 377     pZdauTrue[iDau]=0.;
 378   }
 379
 380   Int_t nDau=partD->GetNDaughters();
 381   Int_t labelFirstDau = partD->GetDaughter(0);
 382   if(nDau==3){
 383     for(Int_t iDau=0; iDau<3; iDau++){
 384       Int_t ind = labelFirstDau+iDau;
 385       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
 386       if(!part){
 387         printf("Daughter particle not found in MC array");
 388         return 99999.;
 389       }
 390       pXdauTrue[iDau]=part->Px();
 391       pYdauTrue[iDau]=part->Py();
 392       pZdauTrue[iDau]=part->Pz();
 393     }
 394   }else if(nDau==2){
 395     Int_t theDau=0;
 396     for(Int_t iDau=0; iDau<2; iDau++){
 397       Int_t ind = labelFirstDau+iDau;
 398       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
 399       if(!part){
 400         printf("Daughter particle not found in MC array");
 401         return 99999.;
 402       }
 403       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
 404       if(pdgCode==211 || pdgCode==321){
 405         pXdauTrue[theDau]=part->Px();
 406         pYdauTrue[theDau]=part->Py();
 407         pZdauTrue[theDau]=part->Pz();
 408         ++theDau;
 409       }else{
 410         Int_t nDauRes=part->GetNDaughters();
 411         if(nDauRes==2){
 412           Int_t labelFirstDauRes = part->GetDaughter(0);
 413           for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
 414             Int_t indDR = labelFirstDauRes+iDauRes;
 415             AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
 416             if(!partDR){
 417               printf("Daughter particle not found in MC array");
 418               return 99999.;
 419             }
 420
 421             Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
 422             if(pdgCodeDR==211 || pdgCodeDR==321){
 423               pXdauTrue[theDau]=partDR->Px();
 424               pYdauTrue[theDau]=partDR->Py();
 425               pZdauTrue[theDau]=partDR->Pz();
 426               ++theDau;
 427             }
 428           }
 429         }
 430       }
 431     }
 432     if(theDau!=3){
 433       printf("Wrong number of decay prongs");
 434       return 99999.;
 435     }
 436   }
 437
 438   Double_t d0dummy[3]={0.,0.,0.};
 439   AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
 440   return aodDvsMC.ImpParXY();
 441
 442 }
 443
 444
 445
 446 //____________________________________________________________________________
 447 Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) {
 448   //
 449   // Correct the number of accepted tracklets based on a TProfile Hist
 450   //
 451   //
 452
 453   if(TMath::Abs(vtxZ)>10.0){
 454     //    printf("ERROR: Z vertex out of range for correction of multiplicity\n");
 455     return uncorrectedNacc;
 456   }
 457
 458   if(!estimatorAvg){
 459     printf("ERROR: Missing TProfile for correction of multiplicity\n");
 460     return uncorrectedNacc;
 461   }
 462
 463   Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
 464
 465   Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
 466
 467   Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
 468
 469   return correctedNacc;
 470 }
 471 //______________________________________________________________________
 472 TString AliVertexingHFUtils::GetGenerator(Int_t label, AliAODMCHeader* header){
 473   // get the name of the generator that produced a given particle
 474
 475   Int_t nsumpart=0;
 476   TList *lh=header->GetCocktailHeaders();
 477   Int_t nh=lh->GetEntries();
 478   for(Int_t i=0;i<nh;i++){
 479     AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
 480     TString genname=gh->GetName();
 481     Int_t npart=gh->NProduced();
 482     if(label>=nsumpart && label<(nsumpart+npart)) return genname;
 483     nsumpart+=npart;
 484   }
 485   TString empty="";
 486   return empty;
 487 }
 488 //_____________________________________________________________________
 489 void AliVertexingHFUtils::GetTrackPrimaryGenerator(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
 490
 491   // method to check if a track comes from a given generator
 492
 493   Int_t lab=TMath::Abs(track->GetLabel());
 494   nameGen=GetGenerator(lab,header);
 495
 496   //  Int_t countControl=0;
 497
 498   while(nameGen.IsWhitespace()){
 499     AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
 500     if(!mcpart){
 501       printf("AliVertexingHFUtils::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
 502       break;
 503     }
 504     Int_t mother = mcpart->GetMother();
 505     if(mother<0){
 506       printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
 507       break;
 508     }
 509     lab=mother;
 510     nameGen=GetGenerator(mother,header);
 511     // countControl++;
 512     // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
 513     //   printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
 514     //   break;
 515     // }
 516   }
 517
 518   return;
 519 }
 520 //----------------------------------------------------------------------
 521 Bool_t AliVertexingHFUtils::IsTrackInjected(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC){
 522   // method to check if a track comes from the signal event or from the underlying Hijing event
 523   TString nameGen;
 524   GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
 525
 526   if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
 527
 528   return kTRUE;
 529 }
 530 //____________________________________________________________________________
 531 Bool_t AliVertexingHFUtils::IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
 532   // method to check if a D meson candidate comes from the signal event or from the underlying Hijing event
 533
 534   Int_t nprongs=cand->GetNProngs();
 535   for(Int_t i=0;i<nprongs;i++){
 536     AliAODTrack *daugh=(AliAODTrack*)cand->GetDaughter(i);
 537     if(IsTrackInjected(daugh,header,arrayMC)) return kTRUE;
 538   }
 539   return kFALSE;
 540 }
 541 //____________________________________________________________________________
 542 Int_t AliVertexingHFUtils::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark){
 543   // checking whether the mother of the particles come from a charm or a bottom quark
 544
 545   Int_t pdgGranma = 0;
 546   Int_t mother = 0;
 547   mother = mcPart->GetMother();
 548   Int_t istep = 0;
 549   Int_t abspdgGranma =0;
 550   Bool_t isFromB=kFALSE;
 551   Bool_t isQuarkFound=kFALSE;
 552   while (mother >0 ){
 553     istep++;
 554     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
 555     if (mcGranma){
 556       pdgGranma = mcGranma->GetPdgCode();
 557       abspdgGranma = TMath::Abs(pdgGranma);
 558       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
 559         isFromB=kTRUE;
 560       }
 561       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
 562       mother = mcGranma->GetMother();
 563     }else{
 564       printf("AliVertexingHFUtils::CheckOrigin: Failed casting the mother particle!");
 565       break;
 566     }
 567   }
 568   if(searchUpToQuark && !isQuarkFound) return 0;
 569   if(isFromB) return 5;
 570   else return 4;
 571
 572 }
 573 //____________________________________________________________________________
 574 Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
 575   // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases
 576
 577
 578   Int_t pdgD=mcPart->GetPdgCode();
 579   if(TMath::Abs(pdgD)!=421) return -1;
 580
 581   Int_t nDau=mcPart->GetNDaughters();
 582
 583   if(nDau==2){
 584     Int_t daughter0 = mcPart->GetDaughter(0);
 585     Int_t daughter1 = mcPart->GetDaughter(1);
 586     AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter0));
 587     AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter1));
 588     if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
 589     arrayDauLab[0]=daughter0;
 590     arrayDauLab[1]=daughter1;
 591     Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
 592     Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
 593     if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
 594        !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
 595       return -1;
 596     }
 597     if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
 598     if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
 599     if((pdgCode0*pdgCode1)>0) return -1;
 600     Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
 601     Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
 602     Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
 603     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
 604     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
 605     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
 606     return 1;
 607   }
 608
 609   if(nDau==3 || nDau==4){
 610     Int_t nKaons=0;
 611     Int_t nPions=0;
 612     Double_t sumPxDau=0.;
 613     Double_t sumPyDau=0.;
 614     Double_t sumPzDau=0.;
 615     Int_t nFoundKpi=0;
 616     Int_t labelFirstDau = mcPart->GetDaughter(0);
 617     for(Int_t iDau=0; iDau<nDau; iDau++){
 618       Int_t indDau = labelFirstDau+iDau;
 619       if(indDau<0) return -1;
 620       AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
 621       if(!dau) return -1;
 622       Int_t pdgdau=dau->GetPdgCode();
 623       if(TMath::Abs(pdgdau)==321){
 624         if(pdgD>0 && pdgdau>0) return -1;
 625         if(pdgD<0 && pdgdau<0) return -1;
 626         nKaons++;
 627         sumPxDau+=dau->Px();
 628         sumPyDau+=dau->Py();
 629         sumPzDau+=dau->Pz();
 630         arrayDauLab[nFoundKpi++]=indDau;
 631         if(nFoundKpi>4) return -1;
 632       }else if(TMath::Abs(pdgdau)==211){
 633         nPions++;
 634         sumPxDau+=dau->Px();
 635         sumPyDau+=dau->Py();
 636         sumPzDau+=dau->Pz();
 637         arrayDauLab[nFoundKpi++]=indDau;
 638         if(nFoundKpi>4) return -1;
 639       }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
 640         Int_t nResDau=dau->GetNDaughters();
 641         if(nResDau!=2) return -1;
 642         Int_t indFirstResDau=dau->GetDaughter(0);
 643         for(Int_t resDau=0; resDau<2; resDau++){
 644           Int_t indResDau=indFirstResDau+resDau;
 645           if(indResDau<0) return -1;
 646           AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
 647           if(!resdau) return -1;
 648           Int_t pdgresdau=resdau->GetPdgCode();
 649           if(TMath::Abs(pdgresdau)==321){
 650             if(pdgD>0 && pdgresdau>0) return -1;
 651             if(pdgD<0 && pdgresdau<0) return -1;
 652             nKaons++;
 653             sumPxDau+=resdau->Px();
 654             sumPyDau+=resdau->Py();
 655             sumPzDau+=resdau->Pz();
 656             arrayDauLab[nFoundKpi++]=indResDau;
 657             if(nFoundKpi>4) return -1;
 658           }
 659           if(TMath::Abs(pdgresdau)==211){
 660             nPions++;
 661             sumPxDau+=resdau->Px();
 662             sumPyDau+=resdau->Py();
 663             sumPzDau+=resdau->Pz();
 664             arrayDauLab[nFoundKpi++]=indResDau;
 665             if(nFoundKpi>4) return -1;
 666           }
 667         }
 668       }else{
 669         return -1;
 670       }
 671     }
 672     if(nPions!=3) return -1;
 673     if(nKaons!=1) return -1;
 674     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
 675     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
 676     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
 677     return 2;
 678   }
 679   return -1;
 680 }
 681 //____________________________________________________________________________
 682 Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
 683   // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases
 684
 685   Int_t pdgD=mcPart->GetPdgCode();
 686   if(TMath::Abs(pdgD)!=411) return -1;
 687
 688   Int_t nDau=mcPart->GetNDaughters();
 689   Int_t labelFirstDau = mcPart->GetDaughter(0);
 690   Int_t nKaons=0;
 691   Int_t nPions=0;
 692   Double_t sumPxDau=0.;
 693   Double_t sumPyDau=0.;
 694   Double_t sumPzDau=0.;
 695   Int_t nFoundKpi=0;
 696
 697   if(nDau==3 || nDau==2){
 698     for(Int_t iDau=0; iDau<nDau; iDau++){
 699       Int_t indDau = labelFirstDau+iDau;
 700       if(indDau<0) return -1;
 701       AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
 702       if(!dau) return -1;
 703       Int_t pdgdau=dau->GetPdgCode();
 704       if(TMath::Abs(pdgdau)==321){
 705         if(pdgD*pdgdau>0) return -1;
 706         nKaons++;
 707         sumPxDau+=dau->Px();
 708         sumPyDau+=dau->Py();
 709         sumPzDau+=dau->Pz();
 710         arrayDauLab[nFoundKpi++]=indDau;
 711         if(nFoundKpi>3) return -1;
 712       }else if(TMath::Abs(pdgdau)==211){
 713         if(pdgD*pdgdau<0) return -1;
 714         nPions++;
 715         sumPxDau+=dau->Px();
 716         sumPyDau+=dau->Py();
 717         sumPzDau+=dau->Pz();
 718         arrayDauLab[nFoundKpi++]=indDau;
 719         if(nFoundKpi>3) return -1;
 720       }else if(TMath::Abs(pdgdau)==313){
 721         Int_t nResDau=dau->GetNDaughters();
 722         if(nResDau!=2) return -1;
 723         Int_t indFirstResDau=dau->GetDaughter(0);
 724         for(Int_t resDau=0; resDau<2; resDau++){
 725           Int_t indResDau=indFirstResDau+resDau;
 726           if(indResDau<0) return -1;
 727           AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
 728           if(!resdau) return -1;
 729           Int_t pdgresdau=resdau->GetPdgCode();
 730           if(TMath::Abs(pdgresdau)==321){
 731             if(pdgD*pdgresdau>0) return -1;
 732             sumPxDau+=resdau->Px();
 733             sumPyDau+=resdau->Py();
 734             sumPzDau+=resdau->Pz();
 735             nKaons++;
 736             arrayDauLab[nFoundKpi++]=indResDau;
 737             if(nFoundKpi>3) return -1;
 738           }
 739           if(TMath::Abs(pdgresdau)==211){
 740             if(pdgD*pdgresdau<0) return -1;
 741             sumPxDau+=resdau->Px();
 742             sumPyDau+=resdau->Py();
 743             sumPzDau+=resdau->Pz();
 744             nPions++;
 745             arrayDauLab[nFoundKpi++]=indResDau;
 746             if(nFoundKpi>3) return -1;
 747           }
 748         }
 749        }else{
 750         return -1;
 751       }
 752     }
 753     if(nPions!=2) return -1;
 754     if(nKaons!=1) return -1;
 755     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
 756     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
 757     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
 758     if(nDau==3) return 1;
 759     else if(nDau==2) return 2;
 760   }
 761
 762   return -1;
 763
 764 }