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 "AliAODMCParticle.h"
  26 #include "AliAODRecoDecayHF.h"
  27 #include "AliVertexingHFUtils.h"
  28
  29 /* $Id$ */
  30
  31 ///////////////////////////////////////////////////////////////////
  32 //                                                               //
  33 // Class with functions useful for different D2H analyses        //
  34 // - event plane resolution                                      //
  35 // - <pt> calculation with side band subtraction                 //
  36 // - tracklet multiplicity calculation                            //
  37 // Origin: F.Prino, Torino, prino@to.infn.it                     //
  38 //                                                               //
  39 ///////////////////////////////////////////////////////////////////
  40
  41 ClassImp(AliVertexingHFUtils)
  42
  43 //______________________________________________________________________
  44 AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
  45   fK(1),
  46   fSubRes(1.),
  47   fMinEtaForTracklets(-1.),
  48   fMaxEtaForTracklets(1.)
  49 {
  50   // Default contructor
  51 }
  52
  53
  54 //______________________________________________________________________
  55 AliVertexingHFUtils::AliVertexingHFUtils(Int_t k):
  56   TObject(),
  57   fK(k),
  58   fSubRes(1.),
  59   fMinEtaForTracklets(-1.),
  60   fMaxEtaForTracklets(1.)
  61 {
  62   // Standard constructor
  63 }
  64
  65
  66 //______________________________________________________________________
  67 Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){
  68   // compute chi from polynomial approximation
  69   Double_t c[5];
  70   if(k==1){
  71     c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
  72   }
  73   else if(k==2){
  74     c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
  75   }
  76   else return -1;
  77   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;
  78 }
  79
  80 //______________________________________________________________________
  81 Double_t AliVertexingHFUtils:: ResolK1(Double_t x){
  82   return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
  83 }
  84
  85
  86 //______________________________________________________________________
  87 Double_t AliVertexingHFUtils::FindChi(Double_t res,  Int_t k){
  88   // compute chi variable (=v2*sqrt(N)) from external values
  89
  90   Double_t x1=0;
  91   Double_t x2=15;
  92   Double_t y1,y2;
  93   if(k==1){
  94     y1=ResolK1(x1)-res;
  95     y2=ResolK1(x2)-res;
  96   }
  97   else if(k==2){
  98     y1=Pol(x1,2)-res;
  99     y2=Pol(x2,2)-res;
 100   }
 101   else return -1;
 102
 103   if(y1*y2>0) return -1;
 104   if(y1==0) return y1;
 105   if(y2==0) return y2;
 106   Double_t xmed,ymed;
 107   Int_t jiter=0;
 108   while((x2-x1)>0.0001){
 109     xmed=0.5*(x1+x2);
 110     if(k==1){
 111       y1=ResolK1(x1)-res;
 112       y2=ResolK1(x2)-res;
 113       ymed=ResolK1(xmed)-res;
 114     }
 115     else if(k==2){
 116       y1=Pol(x1,2)-res;
 117       y2=Pol(x2,2)-res;
 118       ymed=Pol(xmed,2)-res;
 119     }
 120     else return -1;
 121     if((y1*ymed)<0) x2=xmed;
 122     if((y2*ymed)<0)x1=xmed;
 123     if(ymed==0) return xmed;
 124     jiter++;
 125   }
 126   return 0.5*(x1+x2);
 127 }
 128
 129 //______________________________________________________________________
 130 Double_t AliVertexingHFUtils::GetFullEvResol(Double_t resSub, Int_t k){
 131   // computes event plane resolution starting from sub event resolution
 132   Double_t chisub=FindChi(resSub,k);
 133   Double_t chifull=chisub*TMath::Sqrt(2);
 134   if(k==1) return ResolK1(chifull);
 135   else if(k==2) return Pol(chifull,2);
 136   else{
 137     printf("k should be <=2\n");
 138     return 1.;
 139   }
 140 }
 141
 142 //______________________________________________________________________
 143 Double_t AliVertexingHFUtils::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){
 144   // computes event plane resolution starting from sub event correlation histogram
 145   if(!hSubEvCorr) return 1.;
 146   Double_t resSub=GetSubEvResol(hSubEvCorr);
 147   return GetFullEvResol(resSub,k);
 148 }
 149 //______________________________________________________________________
 150 Double_t AliVertexingHFUtils::GetFullEvResolLowLim(const TH1F* hSubEvCorr, Int_t k){
 151   // computes low limit event plane resolution starting from sub event correlation histogram
 152   if(!hSubEvCorr) return 1.;
 153   Double_t resSub=GetSubEvResolLowLim(hSubEvCorr);
 154   printf("%f\n",resSub);
 155   return GetFullEvResol(resSub,k);
 156 }
 157 //______________________________________________________________________
 158 Double_t AliVertexingHFUtils::GetFullEvResolHighLim(const TH1F* hSubEvCorr, Int_t k){
 159   // computes high limit event plane resolution starting from sub event correlation histogram
 160   if(!hSubEvCorr) return 1.;
 161   Double_t resSub=GetSubEvResolHighLim(hSubEvCorr);
 162   printf("%f\n",resSub);
 163   return GetFullEvResol(resSub,k);
 164 }
 165 //______________________________________________________________________
 166 Int_t AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(AliAODEvent* ev, Double_t mineta, Double_t maxeta){
 167   // counts tracklets in given eta range
 168   AliAODTracklets* tracklets=ev->GetTracklets();
 169   Int_t nTr=tracklets->GetNumberOfTracklets();
 170   Int_t count=0;
 171   for(Int_t iTr=0; iTr<nTr; iTr++){
 172     Double_t theta=tracklets->GetTheta(iTr);
 173     Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
 174     if(eta>mineta && eta<maxeta) count++;
 175   }
 176   return count;
 177 }
 178 //______________________________________________________________________
 179 Int_t AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
 180   // counts generated particles in fgiven eta range
 181
 182   Int_t nChargedMC=0;
 183   for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
 184     AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
 185     Int_t charge = part->Charge();
 186     Double_t eta = part->Eta();
 187     if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++;
 188   }
 189   return nChargedMC;
 190 }
 191 //______________________________________________________________________
 192 Int_t AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
 193   // counts generated primary particles in given eta range
 194
 195   Int_t nChargedMC=0;
 196   for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
 197     AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
 198     Int_t charge = part->Charge();
 199     Double_t eta = part->Eta();
 200     if(charge!=0 && eta>mineta && eta<maxeta){
 201       if(part->IsPrimary())nChargedMC++;
 202     }
 203   }
 204   return nChargedMC;
 205 }
 206 //______________________________________________________________________
 207 Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
 208   // counts generated primary particles in given eta range
 209
 210   Int_t nChargedMC=0;
 211   for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
 212     AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
 213     Int_t charge = part->Charge();
 214     Double_t eta = part->Eta();
 215     if(charge!=0 && eta>mineta && eta<maxeta){
 216       if(part->IsPhysicalPrimary())nChargedMC++;
 217     }
 218   }
 219   return nChargedMC;
 220 }
 221 //______________________________________________________________________
 222 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, Int_t rebin){
 223
 224   // Compute <pt> from 2D histogram M vs pt
 225
 226   //Make 2D histos in the desired pt range
 227   Int_t start=hMassD->FindBin(ptmin);
 228   Int_t end=hMassD->FindBin(ptmax)-1;
 229   const Int_t nx=end-start;
 230   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()));
 231   for(Int_t ix=start;ix<end;ix++){
 232     for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){
 233       hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy));
 234       hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy));
 235     }
 236   }
 237
 238   Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
 239   Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
 240   Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig);
 241   Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig);
 242   Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig);
 243   Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig);
 244   //  printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
 245
 246   Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
 247   Int_t minBinBkgLow=2;
 248   Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
 249   Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
 250   Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
 251   Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
 252   Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
 253   Int_t maxBinBkgHi=hMassD->GetNbinsY()-1;
 254   Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
 255   Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
 256   //  printf("BKG Fit Limits = %f %f  && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
 257
 258   Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin);
 259   Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
 260   Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
 261   //  printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
 262
 263   TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow);
 264   TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi);
 265   TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig);
 266
 267   hMptBkgLo->Rebin(rebin);
 268   hMptBkgHi->Rebin(rebin);
 269   hMptSigReg->Rebin(rebin);
 270
 271   hMptBkgLo->Sumw2();
 272   hMptBkgHi->Sumw2();
 273   TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
 274   hMptBkgLoScal->Scale(bkgSig/bkgLow);
 275   TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
 276   hMptBkgHiScal->Scale(bkgSig/bkgHi);
 277
 278   TH1F* hMptBkgAver=0x0;
 279   hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin");
 280   hMptBkgAver->Add(hMptBkgHiScal);
 281   hMptBkgAver->Scale(0.5);
 282   TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin");
 283   hMptSig->Add(hMptBkgAver,-1.);
 284
 285   averagePt = hMptSig->GetMean();
 286   errorPt = hMptSig->GetMeanError();
 287
 288   delete hMptBkgLo;
 289   delete hMptBkgHi;
 290   delete hMptSigReg;
 291   delete hMptBkgLoScal;
 292   delete hMptBkgHiScal;
 293   delete hMptBkgAver;
 294   delete hMassDpt;
 295   delete hMptSig;
 296
 297 }
 298 //____________________________________________________________________________
 299 Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
 300   // true impact parameter calculation for Dzero
 301
 302   if(!partD || !arrayMC || !mcHeader) return 99999.;
 303   Int_t code=partD->GetPdgCode();
 304   if(TMath::Abs(code)!=421) return 99999.;
 305
 306   Double_t vtxTrue[3];
 307   mcHeader->GetVertex(vtxTrue);
 308   Double_t origD[3];
 309   partD->XvYvZv(origD);
 310   Short_t charge=partD->Charge();
 311   Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
 312   for(Int_t iDau=0; iDau<2; iDau++){
 313     pXdauTrue[iDau]=0.;
 314     pYdauTrue[iDau]=0.;
 315     pZdauTrue[iDau]=0.;
 316   }
 317
 318   Int_t nDau=partD->GetNDaughters();
 319   Int_t labelFirstDau = partD->GetDaughter(0);
 320   if(nDau==2){
 321     for(Int_t iDau=0; iDau<2; iDau++){
 322       Int_t ind = labelFirstDau+iDau;
 323       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
 324       if(!part){
 325         printf("Daughter particle not found in MC array");
 326         return 99999.;
 327       }
 328       pXdauTrue[iDau]=part->Px();
 329       pYdauTrue[iDau]=part->Py();
 330       pZdauTrue[iDau]=part->Pz();
 331     }
 332   }else{
 333     return 99999.;
 334   }
 335
 336   Double_t d0dummy[3]={0.,0.,0.};
 337   AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
 338   return aodDvsMC.ImpParXY();
 339
 340 }
 341 //____________________________________________________________________________
 342 Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
 343   // true impact parameter calculation for Dplus
 344
 345   if(!partD || !arrayMC || !mcHeader) return 99999.;
 346   Int_t code=partD->GetPdgCode();
 347   if(TMath::Abs(code)!=411) return 99999.;
 348
 349   Double_t vtxTrue[3];
 350   mcHeader->GetVertex(vtxTrue);
 351   Double_t origD[3];
 352   partD->XvYvZv(origD);
 353   Short_t charge=partD->Charge();
 354   Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
 355   for(Int_t iDau=0; iDau<3; iDau++){
 356     pXdauTrue[iDau]=0.;
 357     pYdauTrue[iDau]=0.;
 358     pZdauTrue[iDau]=0.;
 359   }
 360
 361   Int_t nDau=partD->GetNDaughters();
 362   Int_t labelFirstDau = partD->GetDaughter(0);
 363   if(nDau==3){
 364     for(Int_t iDau=0; iDau<3; iDau++){
 365       Int_t ind = labelFirstDau+iDau;
 366       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
 367       if(!part){
 368         printf("Daughter particle not found in MC array");
 369         return 99999.;
 370       }
 371       pXdauTrue[iDau]=part->Px();
 372       pYdauTrue[iDau]=part->Py();
 373       pZdauTrue[iDau]=part->Pz();
 374     }
 375   }else if(nDau==2){
 376     Int_t theDau=0;
 377     for(Int_t iDau=0; iDau<2; iDau++){
 378       Int_t ind = labelFirstDau+iDau;
 379       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
 380       if(!part){
 381         printf("Daughter particle not found in MC array");
 382         return 99999.;
 383       }
 384       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
 385       if(pdgCode==211 || pdgCode==321){
 386         pXdauTrue[theDau]=part->Px();
 387         pYdauTrue[theDau]=part->Py();
 388         pZdauTrue[theDau]=part->Pz();
 389         ++theDau;
 390       }else{
 391         Int_t nDauRes=part->GetNDaughters();
 392         if(nDauRes==2){
 393           Int_t labelFirstDauRes = part->GetDaughter(0);
 394           for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
 395             Int_t indDR = labelFirstDauRes+iDauRes;
 396             AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
 397             if(!partDR){
 398               printf("Daughter particle not found in MC array");
 399               return 99999.;
 400             }
 401
 402             Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
 403             if(pdgCodeDR==211 || pdgCodeDR==321){
 404               pXdauTrue[theDau]=partDR->Px();
 405               pYdauTrue[theDau]=partDR->Py();
 406               pZdauTrue[theDau]=partDR->Pz();
 407               ++theDau;
 408             }
 409           }
 410         }
 411       }
 412     }
 413     if(theDau!=3){
 414       printf("Wrong number of decay prongs");
 415       return 99999.;
 416     }
 417   }
 418
 419   Double_t d0dummy[3]={0.,0.,0.};
 420   AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
 421   return aodDvsMC.ImpParXY();
 422
 423 }
 424
 425
 426
 427 //____________________________________________________________________________
 428 Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) {
 429   //
 430   // Correct the number of accepted tracklets based on a TProfile Hist
 431   //
 432   //
 433
 434   if(TMath::Abs(vtxZ)>10.0){
 435     printf("ERROR: Z vertex out of range for correction of multiplicity\n");
 436     return uncorrectedNacc;
 437   }
 438
 439   if(!estimatorAvg){
 440     printf("ERROR: Missing TProfile for correction of multiplicity\n");
 441     return uncorrectedNacc;
 442   }
 443
 444   Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
 445
 446   Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
 447
 448   Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
 449
 450   return correctedNacc;
 451 }
 452