Task to evaluate PID efficiency from V0s
[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>
686e4710 17#include <TRandom.h>
a6c5a2e9 18#include "AliAODEvent.h"
25504150 19#include "AliAODMCParticle.h"
a6c5a2e9 20#include "AliVertexingHFUtils.h"
21
d947ea9c 22/* $Id$ */
a6c5a2e9 23
24///////////////////////////////////////////////////////////////////
25// //
26// Class with functions useful for different D2H analyses //
27// - event plane resolution //
28// - <pt> calculation with side band subtraction //
25504150 29// - tracklet multiplicity calculation //
a6c5a2e9 30// Origin: F.Prino, Torino, prino@to.infn.it //
31// //
32///////////////////////////////////////////////////////////////////
33
34ClassImp(AliVertexingHFUtils)
35
36//______________________________________________________________________
37AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
38 fK(1),
39 fSubRes(1.),
40 fMinEtaForTracklets(-1.),
41 fMaxEtaForTracklets(1.)
42{
43 // Default contructor
44}
45
46
47//______________________________________________________________________
48AliVertexingHFUtils::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//______________________________________________________________________
60Double_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//______________________________________________________________________
74Double_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//______________________________________________________________________
80Double_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//______________________________________________________________________
123Double_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//______________________________________________________________________
136Double_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//______________________________________________________________________
143Double_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//______________________________________________________________________
151Double_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//______________________________________________________________________
159Int_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}
86d84fa4 171//______________________________________________________________________
25504150 172Int_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//______________________________________________________________________
185Int_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//______________________________________________________________________
200Int_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//______________________________________________________________________
86d84fa4 215void 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}
686e4710 291//____________________________________________________________________________
292Double_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