1 /**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 #include "AliAODEvent.h"
18 #include "AliAODMCParticle.h"
19 #include "AliVertexingHFUtils.h"
23 ///////////////////////////////////////////////////////////////////
25 // Class with functions useful for different D2H analyses //
26 // - event plane resolution //
27 // - <pt> calculation with side band subtraction //
28 // - tracklet multiplicity calculation //
29 // Origin: F.Prino, Torino, prino@to.infn.it //
31 ///////////////////////////////////////////////////////////////////
33 ClassImp(AliVertexingHFUtils)
35 //______________________________________________________________________
36 AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
39 fMinEtaForTracklets(-1.),
40 fMaxEtaForTracklets(1.)
46 //______________________________________________________________________
47 AliVertexingHFUtils::AliVertexingHFUtils(Int_t k):
51 fMinEtaForTracklets(-1.),
52 fMaxEtaForTracklets(1.)
54 // Standard constructor
58 //______________________________________________________________________
59 Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){
60 // compute chi from polynomial approximation
63 c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
66 c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
69 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;
72 //______________________________________________________________________
73 Double_t AliVertexingHFUtils:: ResolK1(Double_t x){
74 return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
78 //______________________________________________________________________
79 Double_t AliVertexingHFUtils::FindChi(Double_t res, Int_t k){
80 // compute chi variable (=v2*sqrt(N)) from external values
95 if(y1*y2>0) return -1;
100 while((x2-x1)>0.0001){
105 ymed=ResolK1(xmed)-res;
110 ymed=Pol(xmed,2)-res;
113 if((y1*ymed)<0) x2=xmed;
114 if((y2*ymed)<0)x1=xmed;
115 if(ymed==0) return xmed;
121 //______________________________________________________________________
122 Double_t AliVertexingHFUtils::GetFullEvResol(Double_t resSub, Int_t k){
123 // computes event plane resolution starting from sub event resolution
124 Double_t chisub=FindChi(resSub,k);
125 Double_t chifull=chisub*TMath::Sqrt(2);
126 if(k==1) return ResolK1(chifull);
127 else if(k==2) return Pol(chifull,2);
129 printf("k should be <=2\n");
134 //______________________________________________________________________
135 Double_t AliVertexingHFUtils::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){
136 // computes event plane resolution starting from sub event correlation histogram
137 if(!hSubEvCorr) return 1.;
138 Double_t resSub=GetSubEvResol(hSubEvCorr);
139 return GetFullEvResol(resSub,k);
141 //______________________________________________________________________
142 Double_t AliVertexingHFUtils::GetFullEvResolLowLim(const TH1F* hSubEvCorr, Int_t k){
143 // computes low limit event plane resolution starting from sub event correlation histogram
144 if(!hSubEvCorr) return 1.;
145 Double_t resSub=GetSubEvResolLowLim(hSubEvCorr);
146 printf("%f\n",resSub);
147 return GetFullEvResol(resSub,k);
149 //______________________________________________________________________
150 Double_t AliVertexingHFUtils::GetFullEvResolHighLim(const TH1F* hSubEvCorr, Int_t k){
151 // computes high limit event plane resolution starting from sub event correlation histogram
152 if(!hSubEvCorr) return 1.;
153 Double_t resSub=GetSubEvResolHighLim(hSubEvCorr);
154 printf("%f\n",resSub);
155 return GetFullEvResol(resSub,k);
157 //______________________________________________________________________
158 Int_t AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(AliAODEvent* ev, Double_t mineta, Double_t maxeta){
159 // counts tracklets in given eta range
160 AliAODTracklets* tracklets=ev->GetTracklets();
161 Int_t nTr=tracklets->GetNumberOfTracklets();
163 for(Int_t iTr=0; iTr<nTr; iTr++){
164 Double_t theta=tracklets->GetTheta(iTr);
165 Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
166 if(eta>mineta && eta<maxeta) count++;
170 //______________________________________________________________________
171 Int_t AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
172 // counts generated particles in fgiven eta range
175 for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
176 AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
177 Int_t charge = part->Charge();
178 Double_t eta = part->Eta();
179 if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++;
183 //______________________________________________________________________
184 Int_t AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
185 // counts generated primary particles in given eta range
188 for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
189 AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
190 Int_t charge = part->Charge();
191 Double_t eta = part->Eta();
192 if(charge!=0 && eta>mineta && eta<maxeta){
193 if(part->IsPrimary())nChargedMC++;
198 //______________________________________________________________________
199 Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
200 // counts generated primary particles in given eta range
203 for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
204 AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
205 Int_t charge = part->Charge();
206 Double_t eta = part->Eta();
207 if(charge!=0 && eta>mineta && eta<maxeta){
208 if(part->IsPhysicalPrimary())nChargedMC++;
213 //______________________________________________________________________
214 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 // Compute <pt> from 2D histogram M vs pt
218 //Make 2D histos in the desired pt range
219 Int_t start=hMassD->FindBin(ptmin);
220 Int_t end=hMassD->FindBin(ptmax)-1;
221 const Int_t nx=end-start;
222 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()));
223 for(Int_t ix=start;ix<end;ix++){
224 for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){
225 hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy));
226 hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy));
230 Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
231 Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
232 Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig);
233 Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig);
234 Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig);
235 Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig);
236 // printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
238 Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
239 Int_t minBinBkgLow=2;
240 Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
241 Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
242 Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
243 Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
244 Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
245 Int_t maxBinBkgHi=hMassD->GetNbinsY()-1;
246 Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
247 Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
248 // printf("BKG Fit Limits = %f %f && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
250 Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin);
251 Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
252 Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
253 // printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
255 TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow);
256 TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi);
257 TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig);
259 hMptBkgLo->Rebin(rebin);
260 hMptBkgHi->Rebin(rebin);
261 hMptSigReg->Rebin(rebin);
265 TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
266 hMptBkgLoScal->Scale(bkgSig/bkgLow);
267 TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
268 hMptBkgHiScal->Scale(bkgSig/bkgHi);
270 TH1F* hMptBkgAver=0x0;
271 hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin");
272 hMptBkgAver->Add(hMptBkgHiScal);
273 hMptBkgAver->Scale(0.5);
274 TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin");
275 hMptSig->Add(hMptBkgAver,-1.);
277 averagePt = hMptSig->GetMean();
278 errorPt = hMptSig->GetMeanError();
283 delete hMptBkgLoScal;
284 delete hMptBkgHiScal;