PWG3/vertexingHF/charmFlow/AliEventPlaneResolution.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 "AliEventPlaneResolution.h"
  18
  19 /* $Id$ */
  20
  21 ///////////////////////////////////////////////////////////////////
  22 //                                                               //
  23 // Implementation of the class to compute event plane resolution //
  24 // for flow analyses                                             //
  25 // Origin: F.Prino, Torino, prino@to.infn.it                     //
  26 //                                                               //
  27 ///////////////////////////////////////////////////////////////////
  28
  29 ClassImp(AliEventPlaneResolution)
  30
  31 //______________________________________________________________________
  32 AliEventPlaneResolution::AliEventPlaneResolution():TObject(),
  33   fK(1),
  34   fSubRes(1.)
  35 {
  36   // Default contructor
  37 }
  38
  39
  40 //______________________________________________________________________
  41 AliEventPlaneResolution::AliEventPlaneResolution(Int_t k):
  42   TObject(),
  43   fK(k),
  44   fSubRes(1.)
  45 {
  46   // Standard constructor
  47 }
  48
  49
  50 //______________________________________________________________________
  51 Double_t AliEventPlaneResolution::Pol(Double_t x, Int_t k){
  52   // compute chi from polynomial approximation
  53   Double_t c[5];
  54   if(k==1){
  55     c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
  56   }
  57   else if(k==2){
  58     c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
  59   }
  60   else return -1;
  61   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;
  62 }
  63
  64 //______________________________________________________________________
  65 Double_t AliEventPlaneResolution:: ResolK1(Double_t x){
  66   return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
  67 }
  68
  69
  70 //______________________________________________________________________
  71 Double_t AliEventPlaneResolution::FindChi(Double_t res,  Int_t k){
  72   // compute chi variable (=v2*sqrt(N)) from external values
  73
  74   Double_t x1=0;
  75   Double_t x2=15;
  76   Double_t y1,y2;
  77   if(k==1){
  78     y1=ResolK1(x1)-res;
  79     y2=ResolK1(x2)-res;
  80   }
  81   else if(k==2){
  82     y1=Pol(x1,2)-res;
  83     y2=Pol(x2,2)-res;
  84   }
  85   else return -1;
  86
  87   if(y1*y2>0) return -1;
  88   if(y1==0) return y1;
  89   if(y2==0) return y2;
  90   Double_t xmed,ymed;
  91   Int_t jiter=0;
  92   while((x2-x1)>0.0001){
  93     xmed=0.5*(x1+x2);
  94     if(k==1){
  95       y1=ResolK1(x1)-res;
  96       y2=ResolK1(x2)-res;
  97       ymed=ResolK1(xmed)-res;
  98     }
  99     else if(k==2){
 100       y1=Pol(x1,2)-res;
 101       y2=Pol(x2,2)-res;
 102       ymed=Pol(xmed,2)-res;
 103     }
 104     else return -1;
 105     if((y1*ymed)<0) x2=xmed;
 106     if((y2*ymed)<0)x1=xmed;
 107     if(ymed==0) return xmed;
 108     jiter++;
 109   }
 110   return 0.5*(x1+x2);
 111 }
 112
 113 //______________________________________________________________________
 114 Double_t AliEventPlaneResolution::GetFullEvResol(Double_t resSub, Int_t k){
 115   // computes event plane resolution starting from sub event resolution
 116   Double_t chisub=FindChi(resSub,k);
 117   Double_t chifull=chisub*TMath::Sqrt(2);
 118   if(k==1) return ResolK1(chifull);
 119   else if(k==2) return Pol(chifull,2);
 120   else{
 121     printf("k should be <=2\n");
 122     return 1.;
 123   }
 124 }
 125
 126 //______________________________________________________________________
 127 Double_t AliEventPlaneResolution::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){
 128   // computes event plane resolution starting from sub event correlation histogram
 129   if(!hSubEvCorr) return 1.;
 130   Double_t resSub=GetSubEvResol(hSubEvCorr);
 131   return GetFullEvResol(resSub,k);
 132 }