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 "AliAODEvent.h"
  18 #include "AliVertexingHFUtils.h"
  19
  20 /* $Id: $ */
  21
  22 ///////////////////////////////////////////////////////////////////
  23 //                                                               //
  24 // Class with functions useful for different D2H analyses        //
  25 // - event plane resolution                                      //
  26 // - <pt> calculation with side band subtraction                 //
  27 // - tracklet multiplcity calculation                            //
  28 // Origin: F.Prino, Torino, prino@to.infn.it                     //
  29 //                                                               //
  30 ///////////////////////////////////////////////////////////////////
  31
  32 ClassImp(AliVertexingHFUtils)
  33
  34 //______________________________________________________________________
  35 AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
  36   fK(1),
  37   fSubRes(1.),
  38   fMinEtaForTracklets(-1.),
  39   fMaxEtaForTracklets(1.)
  40 {
  41   // Default contructor
  42 }
  43
  44
  45 //______________________________________________________________________
  46 AliVertexingHFUtils::AliVertexingHFUtils(Int_t k):
  47   TObject(),
  48   fK(k),
  49   fSubRes(1.),
  50   fMinEtaForTracklets(-1.),
  51   fMaxEtaForTracklets(1.)
  52 {
  53   // Standard constructor
  54 }
  55
  56
  57 //______________________________________________________________________
  58 Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){
  59   // compute chi from polynomial approximation
  60   Double_t c[5];
  61   if(k==1){
  62     c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
  63   }
  64   else if(k==2){
  65     c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
  66   }
  67   else return -1;
  68   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;
  69 }
  70
  71 //______________________________________________________________________
  72 Double_t AliVertexingHFUtils:: ResolK1(Double_t x){
  73   return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
  74 }
  75
  76
  77 //______________________________________________________________________
  78 Double_t AliVertexingHFUtils::FindChi(Double_t res,  Int_t k){
  79   // compute chi variable (=v2*sqrt(N)) from external values
  80
  81   Double_t x1=0;
  82   Double_t x2=15;
  83   Double_t y1,y2;
  84   if(k==1){
  85     y1=ResolK1(x1)-res;
  86     y2=ResolK1(x2)-res;
  87   }
  88   else if(k==2){
  89     y1=Pol(x1,2)-res;
  90     y2=Pol(x2,2)-res;
  91   }
  92   else return -1;
  93
  94   if(y1*y2>0) return -1;
  95   if(y1==0) return y1;
  96   if(y2==0) return y2;
  97   Double_t xmed,ymed;
  98   Int_t jiter=0;
  99   while((x2-x1)>0.0001){
 100     xmed=0.5*(x1+x2);
 101     if(k==1){
 102       y1=ResolK1(x1)-res;
 103       y2=ResolK1(x2)-res;
 104       ymed=ResolK1(xmed)-res;
 105     }
 106     else if(k==2){
 107       y1=Pol(x1,2)-res;
 108       y2=Pol(x2,2)-res;
 109       ymed=Pol(xmed,2)-res;
 110     }
 111     else return -1;
 112     if((y1*ymed)<0) x2=xmed;
 113     if((y2*ymed)<0)x1=xmed;
 114     if(ymed==0) return xmed;
 115     jiter++;
 116   }
 117   return 0.5*(x1+x2);
 118 }
 119
 120 //______________________________________________________________________
 121 Double_t AliVertexingHFUtils::GetFullEvResol(Double_t resSub, Int_t k){
 122   // computes event plane resolution starting from sub event resolution
 123   Double_t chisub=FindChi(resSub,k);
 124   Double_t chifull=chisub*TMath::Sqrt(2);
 125   if(k==1) return ResolK1(chifull);
 126   else if(k==2) return Pol(chifull,2);
 127   else{
 128     printf("k should be <=2\n");
 129     return 1.;
 130   }
 131 }
 132
 133 //______________________________________________________________________
 134 Double_t AliVertexingHFUtils::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){
 135   // computes event plane resolution starting from sub event correlation histogram
 136   if(!hSubEvCorr) return 1.;
 137   Double_t resSub=GetSubEvResol(hSubEvCorr);
 138   return GetFullEvResol(resSub,k);
 139 }
 140 //______________________________________________________________________
 141 Double_t AliVertexingHFUtils::GetFullEvResolLowLim(const TH1F* hSubEvCorr, Int_t k){
 142   // computes low limit event plane resolution starting from sub event correlation histogram
 143   if(!hSubEvCorr) return 1.;
 144   Double_t resSub=GetSubEvResolLowLim(hSubEvCorr);
 145   printf("%f\n",resSub);
 146   return GetFullEvResol(resSub,k);
 147 }
 148 //______________________________________________________________________
 149 Double_t AliVertexingHFUtils::GetFullEvResolHighLim(const TH1F* hSubEvCorr, Int_t k){
 150   // computes high limit event plane resolution starting from sub event correlation histogram
 151   if(!hSubEvCorr) return 1.;
 152   Double_t resSub=GetSubEvResolHighLim(hSubEvCorr);
 153   printf("%f\n",resSub);
 154   return GetFullEvResol(resSub,k);
 155 }
 156 //______________________________________________________________________
 157 Int_t AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(AliAODEvent* ev, Double_t mineta, Double_t maxeta){
 158   // counts tracklets in given eta range
 159   AliAODTracklets* tracklets=ev->GetTracklets();
 160   Int_t nTr=tracklets->GetNumberOfTracklets();
 161   Int_t count=0;
 162   for(Int_t iTr=0; iTr<nTr; iTr++){
 163     Double_t theta=tracklets->GetTheta(iTr);
 164     Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
 165     if(eta>mineta && eta<maxeta) count++;
 166   }
 167   return count;
 168 }