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 "AliAODEvent.h"
  19 #include "AliAODMCParticle.h"
  20 #include "AliVertexingHFUtils.h"
  21
  22 /* $Id$ */
  23
  24 ///////////////////////////////////////////////////////////////////
  25 //                                                               //
  26 // Class with functions useful for different D2H analyses        //
  27 // - event plane resolution                                      //
  28 // - <pt> calculation with side band subtraction                 //
  29 // - tracklet multiplicity calculation                            //
  30 // Origin: F.Prino, Torino, prino@to.infn.it                     //
  31 //                                                               //
  32 ///////////////////////////////////////////////////////////////////
  33
  34 ClassImp(AliVertexingHFUtils)
  35
  36 //______________________________________________________________________
  37 AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
  38   fK(1),
  39   fSubRes(1.),
  40   fMinEtaForTracklets(-1.),
  41   fMaxEtaForTracklets(1.)
  42 {
  43   // Default contructor
  44 }
  45
  46
  47 //______________________________________________________________________
  48 AliVertexingHFUtils::AliVertexingHFUtils(Int_t k):
  49   TObject(),
  50   fK(k),
  51   fSubRes(1.),
  52   fMinEtaForTracklets(-1.),
  53   fMaxEtaForTracklets(1.)
  54 {
  55   // Standard constructor
  56 }
  57
  58
  59 //______________________________________________________________________
  60 Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){
  61   // compute chi from polynomial approximation
  62   Double_t c[5];
  63   if(k==1){
  64     c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
  65   }
  66   else if(k==2){
  67     c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
  68   }
  69   else return -1;
  70   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;
  71 }
  72
  73 //______________________________________________________________________
  74 Double_t AliVertexingHFUtils:: ResolK1(Double_t x){
  75   return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
  76 }
  77
  78
  79 //______________________________________________________________________
  80 Double_t AliVertexingHFUtils::FindChi(Double_t res,  Int_t k){
  81   // compute chi variable (=v2*sqrt(N)) from external values
  82
  83   Double_t x1=0;
  84   Double_t x2=15;
  85   Double_t y1,y2;
  86   if(k==1){
  87     y1=ResolK1(x1)-res;
  88     y2=ResolK1(x2)-res;
  89   }
  90   else if(k==2){
  91     y1=Pol(x1,2)-res;
  92     y2=Pol(x2,2)-res;
  93   }
  94   else return -1;
  95
  96   if(y1*y2>0) return -1;
  97   if(y1==0) return y1;
  98   if(y2==0) return y2;
  99   Double_t xmed,ymed;
 100   Int_t jiter=0;
 101   while((x2-x1)>0.0001){
 102     xmed=0.5*(x1+x2);
 103     if(k==1){
 104       y1=ResolK1(x1)-res;
 105       y2=ResolK1(x2)-res;
 106       ymed=ResolK1(xmed)-res;
 107     }
 108     else if(k==2){
 109       y1=Pol(x1,2)-res;
 110       y2=Pol(x2,2)-res;
 111       ymed=Pol(xmed,2)-res;
 112     }
 113     else return -1;
 114     if((y1*ymed)<0) x2=xmed;
 115     if((y2*ymed)<0)x1=xmed;
 116     if(ymed==0) return xmed;
 117     jiter++;
 118   }
 119   return 0.5*(x1+x2);
 120 }
 121
 122 //______________________________________________________________________
 123 Double_t AliVertexingHFUtils::GetFullEvResol(Double_t resSub, Int_t k){
 124   // computes event plane resolution starting from sub event resolution
 125   Double_t chisub=FindChi(resSub,k);
 126   Double_t chifull=chisub*TMath::Sqrt(2);
 127   if(k==1) return ResolK1(chifull);
 128   else if(k==2) return Pol(chifull,2);
 129   else{
 130     printf("k should be <=2\n");
 131     return 1.;
 132   }
 133 }
 134
 135 //______________________________________________________________________
 136 Double_t AliVertexingHFUtils::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){
 137   // computes event plane resolution starting from sub event correlation histogram
 138   if(!hSubEvCorr) return 1.;
 139   Double_t resSub=GetSubEvResol(hSubEvCorr);
 140   return GetFullEvResol(resSub,k);
 141 }
 142 //______________________________________________________________________
 143 Double_t AliVertexingHFUtils::GetFullEvResolLowLim(const TH1F* hSubEvCorr, Int_t k){
 144   // computes low limit event plane resolution starting from sub event correlation histogram
 145   if(!hSubEvCorr) return 1.;
 146   Double_t resSub=GetSubEvResolLowLim(hSubEvCorr);
 147   printf("%f\n",resSub);
 148   return GetFullEvResol(resSub,k);
 149 }
 150 //______________________________________________________________________
 151 Double_t AliVertexingHFUtils::GetFullEvResolHighLim(const TH1F* hSubEvCorr, Int_t k){
 152   // computes high limit event plane resolution starting from sub event correlation histogram
 153   if(!hSubEvCorr) return 1.;
 154   Double_t resSub=GetSubEvResolHighLim(hSubEvCorr);
 155   printf("%f\n",resSub);
 156   return GetFullEvResol(resSub,k);
 157 }
 158 //______________________________________________________________________
 159 Int_t AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(AliAODEvent* ev, Double_t mineta, Double_t maxeta){
 160   // counts tracklets in given eta range
 161   AliAODTracklets* tracklets=ev->GetTracklets();
 162   Int_t nTr=tracklets->GetNumberOfTracklets();
 163   Int_t count=0;
 164   for(Int_t iTr=0; iTr<nTr; iTr++){
 165     Double_t theta=tracklets->GetTheta(iTr);
 166     Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
 167     if(eta>mineta && eta<maxeta) count++;
 168   }
 169   return count;
 170 }
 171 //______________________________________________________________________
 172 Int_t AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
 173   // counts generated particles in fgiven eta range
 174
 175   Int_t nChargedMC=0;
 176   for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
 177     AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
 178     Int_t charge = part->Charge();
 179     Double_t eta = part->Eta();
 180     if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++;
 181   }
 182   return nChargedMC;
 183 }
 184 //______________________________________________________________________
 185 Int_t AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
 186   // counts generated primary particles in given eta range
 187
 188   Int_t nChargedMC=0;
 189   for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
 190     AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
 191     Int_t charge = part->Charge();
 192     Double_t eta = part->Eta();
 193     if(charge!=0 && eta>mineta && eta<maxeta){
 194       if(part->IsPrimary())nChargedMC++;
 195     }
 196   }
 197   return nChargedMC;
 198 }
 199 //______________________________________________________________________
 200 Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
 201   // counts generated primary particles in given 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){
 209       if(part->IsPhysicalPrimary())nChargedMC++;
 210     }
 211   }
 212   return nChargedMC;
 213 }
 214 //______________________________________________________________________
 215 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){
 216
 217   // Compute <pt> from 2D histogram M vs pt
 218
 219   //Make 2D histos in the desired pt range
 220   Int_t start=hMassD->FindBin(ptmin);
 221   Int_t end=hMassD->FindBin(ptmax)-1;
 222   const Int_t nx=end-start;
 223   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()));
 224   for(Int_t ix=start;ix<end;ix++){
 225     for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){
 226       hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy));
 227       hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy));
 228     }
 229   }
 230
 231   Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
 232   Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
 233   Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig);
 234   Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig);
 235   Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig);
 236   Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig);
 237   //  printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
 238
 239   Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
 240   Int_t minBinBkgLow=2;
 241   Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
 242   Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
 243   Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
 244   Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
 245   Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
 246   Int_t maxBinBkgHi=hMassD->GetNbinsY()-1;
 247   Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
 248   Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
 249   //  printf("BKG Fit Limits = %f %f  && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
 250
 251   Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin);
 252   Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
 253   Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
 254   //  printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
 255
 256   TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow);
 257   TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi);
 258   TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig);
 259
 260   hMptBkgLo->Rebin(rebin);
 261   hMptBkgHi->Rebin(rebin);
 262   hMptSigReg->Rebin(rebin);
 263
 264   hMptBkgLo->Sumw2();
 265   hMptBkgHi->Sumw2();
 266   TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
 267   hMptBkgLoScal->Scale(bkgSig/bkgLow);
 268   TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
 269   hMptBkgHiScal->Scale(bkgSig/bkgHi);
 270
 271   TH1F* hMptBkgAver=0x0;
 272   hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin");
 273   hMptBkgAver->Add(hMptBkgHiScal);
 274   hMptBkgAver->Scale(0.5);
 275   TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin");
 276   hMptSig->Add(hMptBkgAver,-1.);
 277
 278   averagePt = hMptSig->GetMean();
 279   errorPt = hMptSig->GetMeanError();
 280
 281   delete hMptBkgLo;
 282   delete hMptBkgHi;
 283   delete hMptSigReg;
 284   delete hMptBkgLoScal;
 285   delete hMptBkgHiScal;
 286   delete hMptBkgAver;
 287   delete hMassDpt;
 288   delete hMptSig;
 289
 290 }
 291 //____________________________________________________________________________
 292 Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) {
 293   //
 294   // Correct the number of accepted tracklets based on a TProfile Hist
 295   //
 296   //
 297
 298   if(TMath::Abs(vtxZ)>10.0){
 299     printf("ERROR: Z vertex out of range for correction of multiplicity\n");
 300     return uncorrectedNacc;
 301   }
 302
 303   if(!estimatorAvg){
 304     printf("ERROR: Missing TProfile for correction of multiplicity\n");
 305     return uncorrectedNacc;
 306   }
 307
 308   Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
 309
 310   Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
 311
 312   Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
 313
 314   return correctedNacc;
 315 }
 316