New class for common functions of D2H analyses (not yet included in compilation)
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliVertexingHFUtils.cxx
CommitLineData
a6c5a2e9 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
32ClassImp(AliVertexingHFUtils)
33
34//______________________________________________________________________
35AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
36 fK(1),
37 fSubRes(1.),
38 fMinEtaForTracklets(-1.),
39 fMaxEtaForTracklets(1.)
40{
41 // Default contructor
42}
43
44
45//______________________________________________________________________
46AliVertexingHFUtils::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//______________________________________________________________________
58Double_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//______________________________________________________________________
72Double_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//______________________________________________________________________
78Double_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//______________________________________________________________________
121Double_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//______________________________________________________________________
134Double_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//______________________________________________________________________
141Double_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//______________________________________________________________________
149Double_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//______________________________________________________________________
157Int_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}