]>
Commit | Line | Data |
---|---|---|
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> |
b176db3b | 18 | #include <TProfile.h> |
19 | #include <TClonesArray.h> | |
20 | #include <TH1F.h> | |
21 | #include <TH2F.h> | |
22 | #include <TF1.h> | |
a6c5a2e9 | 23 | #include "AliAODEvent.h" |
b176db3b | 24 | #include "AliAODMCHeader.h" |
25504150 | 25 | #include "AliAODMCParticle.h" |
b176db3b | 26 | #include "AliAODRecoDecayHF.h" |
a6c5a2e9 | 27 | #include "AliVertexingHFUtils.h" |
28 | ||
d947ea9c | 29 | /* $Id$ */ |
a6c5a2e9 | 30 | |
31 | /////////////////////////////////////////////////////////////////// | |
32 | // // | |
33 | // Class with functions useful for different D2H analyses // | |
34 | // - event plane resolution // | |
35 | // - <pt> calculation with side band subtraction // | |
25504150 | 36 | // - tracklet multiplicity calculation // |
a6c5a2e9 | 37 | // Origin: F.Prino, Torino, prino@to.infn.it // |
38 | // // | |
39 | /////////////////////////////////////////////////////////////////// | |
40 | ||
41 | ClassImp(AliVertexingHFUtils) | |
42 | ||
43 | //______________________________________________________________________ | |
44 | AliVertexingHFUtils::AliVertexingHFUtils():TObject(), | |
45 | fK(1), | |
46 | fSubRes(1.), | |
47 | fMinEtaForTracklets(-1.), | |
48 | fMaxEtaForTracklets(1.) | |
49 | { | |
50 | // Default contructor | |
51 | } | |
52 | ||
53 | ||
54 | //______________________________________________________________________ | |
55 | AliVertexingHFUtils::AliVertexingHFUtils(Int_t k): | |
56 | TObject(), | |
57 | fK(k), | |
58 | fSubRes(1.), | |
59 | fMinEtaForTracklets(-1.), | |
60 | fMaxEtaForTracklets(1.) | |
61 | { | |
62 | // Standard constructor | |
63 | } | |
64 | ||
65 | ||
5d1396f9 | 66 | //______________________________________________________________________ |
67 | void AliVertexingHFUtils::ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance,Double_t &errsignificance){ | |
68 | // calculate significance from S, B and errors | |
69 | ||
70 | ||
71 | Double_t errSigSq=errsignal*errsignal; | |
72 | Double_t errBkgSq=errbackground*errbackground; | |
73 | Double_t sigPlusBkg=signal+background; | |
74 | if(sigPlusBkg>0. && signal>0.){ | |
75 | significance = signal/TMath::Sqrt(signal+background); | |
76 | errsignificance = significance*TMath::Sqrt((errSigSq+errBkgSq)/(4.*sigPlusBkg*sigPlusBkg)+(background/sigPlusBkg)*errSigSq/signal/signal); | |
77 | }else{ | |
78 | significance=0.; | |
79 | errsignificance=0.; | |
80 | } | |
81 | return; | |
82 | ||
83 | } | |
a6c5a2e9 | 84 | //______________________________________________________________________ |
85 | Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){ | |
86 | // compute chi from polynomial approximation | |
87 | Double_t c[5]; | |
88 | if(k==1){ | |
89 | c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283; | |
90 | } | |
91 | else if(k==2){ | |
92 | c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815; | |
93 | } | |
94 | else return -1; | |
95 | 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; | |
96 | } | |
97 | ||
98 | //______________________________________________________________________ | |
99 | Double_t AliVertexingHFUtils:: ResolK1(Double_t x){ | |
100 | return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4)); | |
101 | } | |
102 | ||
103 | ||
104 | //______________________________________________________________________ | |
105 | Double_t AliVertexingHFUtils::FindChi(Double_t res, Int_t k){ | |
106 | // compute chi variable (=v2*sqrt(N)) from external values | |
107 | ||
108 | Double_t x1=0; | |
109 | Double_t x2=15; | |
110 | Double_t y1,y2; | |
111 | if(k==1){ | |
112 | y1=ResolK1(x1)-res; | |
113 | y2=ResolK1(x2)-res; | |
114 | } | |
115 | else if(k==2){ | |
116 | y1=Pol(x1,2)-res; | |
117 | y2=Pol(x2,2)-res; | |
118 | } | |
119 | else return -1; | |
120 | ||
121 | if(y1*y2>0) return -1; | |
122 | if(y1==0) return y1; | |
123 | if(y2==0) return y2; | |
124 | Double_t xmed,ymed; | |
125 | Int_t jiter=0; | |
126 | while((x2-x1)>0.0001){ | |
127 | xmed=0.5*(x1+x2); | |
128 | if(k==1){ | |
129 | y1=ResolK1(x1)-res; | |
130 | y2=ResolK1(x2)-res; | |
131 | ymed=ResolK1(xmed)-res; | |
132 | } | |
133 | else if(k==2){ | |
134 | y1=Pol(x1,2)-res; | |
135 | y2=Pol(x2,2)-res; | |
136 | ymed=Pol(xmed,2)-res; | |
137 | } | |
138 | else return -1; | |
139 | if((y1*ymed)<0) x2=xmed; | |
140 | if((y2*ymed)<0)x1=xmed; | |
141 | if(ymed==0) return xmed; | |
142 | jiter++; | |
143 | } | |
144 | return 0.5*(x1+x2); | |
145 | } | |
146 | ||
147 | //______________________________________________________________________ | |
148 | Double_t AliVertexingHFUtils::GetFullEvResol(Double_t resSub, Int_t k){ | |
149 | // computes event plane resolution starting from sub event resolution | |
150 | Double_t chisub=FindChi(resSub,k); | |
151 | Double_t chifull=chisub*TMath::Sqrt(2); | |
152 | if(k==1) return ResolK1(chifull); | |
153 | else if(k==2) return Pol(chifull,2); | |
154 | else{ | |
155 | printf("k should be <=2\n"); | |
156 | return 1.; | |
157 | } | |
158 | } | |
159 | ||
160 | //______________________________________________________________________ | |
161 | Double_t AliVertexingHFUtils::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){ | |
162 | // computes event plane resolution starting from sub event correlation histogram | |
163 | if(!hSubEvCorr) return 1.; | |
164 | Double_t resSub=GetSubEvResol(hSubEvCorr); | |
165 | return GetFullEvResol(resSub,k); | |
166 | } | |
167 | //______________________________________________________________________ | |
168 | Double_t AliVertexingHFUtils::GetFullEvResolLowLim(const TH1F* hSubEvCorr, Int_t k){ | |
169 | // computes low limit event plane resolution starting from sub event correlation histogram | |
170 | if(!hSubEvCorr) return 1.; | |
171 | Double_t resSub=GetSubEvResolLowLim(hSubEvCorr); | |
172 | printf("%f\n",resSub); | |
173 | return GetFullEvResol(resSub,k); | |
174 | } | |
175 | //______________________________________________________________________ | |
176 | Double_t AliVertexingHFUtils::GetFullEvResolHighLim(const TH1F* hSubEvCorr, Int_t k){ | |
177 | // computes high limit event plane resolution starting from sub event correlation histogram | |
178 | if(!hSubEvCorr) return 1.; | |
179 | Double_t resSub=GetSubEvResolHighLim(hSubEvCorr); | |
180 | printf("%f\n",resSub); | |
181 | return GetFullEvResol(resSub,k); | |
182 | } | |
183 | //______________________________________________________________________ | |
184 | Int_t AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(AliAODEvent* ev, Double_t mineta, Double_t maxeta){ | |
185 | // counts tracklets in given eta range | |
186 | AliAODTracklets* tracklets=ev->GetTracklets(); | |
187 | Int_t nTr=tracklets->GetNumberOfTracklets(); | |
188 | Int_t count=0; | |
189 | for(Int_t iTr=0; iTr<nTr; iTr++){ | |
190 | Double_t theta=tracklets->GetTheta(iTr); | |
191 | Double_t eta=-TMath::Log(TMath::Tan(theta/2.)); | |
192 | if(eta>mineta && eta<maxeta) count++; | |
193 | } | |
194 | return count; | |
195 | } | |
86d84fa4 | 196 | //______________________________________________________________________ |
25504150 | 197 | Int_t AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){ |
198 | // counts generated particles in fgiven eta range | |
199 | ||
200 | Int_t nChargedMC=0; | |
201 | for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){ | |
202 | AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i); | |
203 | Int_t charge = part->Charge(); | |
204 | Double_t eta = part->Eta(); | |
205 | if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++; | |
206 | } | |
207 | return nChargedMC; | |
208 | } | |
209 | //______________________________________________________________________ | |
210 | Int_t AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){ | |
211 | // counts generated primary particles in given eta range | |
212 | ||
213 | Int_t nChargedMC=0; | |
214 | for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){ | |
215 | AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i); | |
216 | Int_t charge = part->Charge(); | |
217 | Double_t eta = part->Eta(); | |
218 | if(charge!=0 && eta>mineta && eta<maxeta){ | |
219 | if(part->IsPrimary())nChargedMC++; | |
220 | } | |
221 | } | |
222 | return nChargedMC; | |
223 | } | |
224 | //______________________________________________________________________ | |
225 | Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){ | |
226 | // counts generated primary particles in given eta range | |
227 | ||
228 | Int_t nChargedMC=0; | |
229 | for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){ | |
230 | AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i); | |
231 | Int_t charge = part->Charge(); | |
232 | Double_t eta = part->Eta(); | |
233 | if(charge!=0 && eta>mineta && eta<maxeta){ | |
234 | if(part->IsPhysicalPrimary())nChargedMC++; | |
235 | } | |
236 | } | |
237 | return nChargedMC; | |
238 | } | |
239 | //______________________________________________________________________ | |
86d84fa4 | 240 | 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){ |
241 | ||
242 | // Compute <pt> from 2D histogram M vs pt | |
243 | ||
244 | //Make 2D histos in the desired pt range | |
245 | Int_t start=hMassD->FindBin(ptmin); | |
246 | Int_t end=hMassD->FindBin(ptmax)-1; | |
247 | const Int_t nx=end-start; | |
248 | 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())); | |
249 | for(Int_t ix=start;ix<end;ix++){ | |
250 | for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){ | |
251 | hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy)); | |
252 | hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy)); | |
253 | } | |
254 | } | |
255 | ||
256 | Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit; | |
257 | Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit; | |
258 | Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig); | |
259 | Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig); | |
260 | Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig); | |
261 | Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig); | |
262 | // printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin); | |
263 | ||
264 | Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit; | |
265 | Int_t minBinBkgLow=2; | |
266 | Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow); | |
267 | Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow); | |
268 | Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow); | |
269 | Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit; | |
270 | Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi); | |
271 | Int_t maxBinBkgHi=hMassD->GetNbinsY()-1; | |
272 | Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi); | |
273 | Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi); | |
274 | // printf("BKG Fit Limits = %f %f && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin); | |
275 | ||
276 | Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin); | |
277 | Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin); | |
278 | Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin); | |
279 | // printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi); | |
280 | ||
281 | TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow); | |
282 | TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi); | |
283 | TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig); | |
284 | ||
285 | hMptBkgLo->Rebin(rebin); | |
286 | hMptBkgHi->Rebin(rebin); | |
287 | hMptSigReg->Rebin(rebin); | |
288 | ||
289 | hMptBkgLo->Sumw2(); | |
290 | hMptBkgHi->Sumw2(); | |
291 | TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin"); | |
292 | hMptBkgLoScal->Scale(bkgSig/bkgLow); | |
293 | TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin"); | |
294 | hMptBkgHiScal->Scale(bkgSig/bkgHi); | |
295 | ||
296 | TH1F* hMptBkgAver=0x0; | |
297 | hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin"); | |
298 | hMptBkgAver->Add(hMptBkgHiScal); | |
299 | hMptBkgAver->Scale(0.5); | |
300 | TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin"); | |
301 | hMptSig->Add(hMptBkgAver,-1.); | |
302 | ||
303 | averagePt = hMptSig->GetMean(); | |
304 | errorPt = hMptSig->GetMeanError(); | |
305 | ||
306 | delete hMptBkgLo; | |
307 | delete hMptBkgHi; | |
308 | delete hMptSigReg; | |
309 | delete hMptBkgLoScal; | |
310 | delete hMptBkgHiScal; | |
311 | delete hMptBkgAver; | |
312 | delete hMassDpt; | |
313 | delete hMptSig; | |
314 | ||
15917f37 | 315 | } |
316 | //____________________________________________________________________________ | |
317 | Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) { | |
318 | // true impact parameter calculation for Dzero | |
319 | ||
320 | if(!partD || !arrayMC || !mcHeader) return 99999.; | |
321 | Int_t code=partD->GetPdgCode(); | |
322 | if(TMath::Abs(code)!=421) return 99999.; | |
323 | ||
324 | Double_t vtxTrue[3]; | |
325 | mcHeader->GetVertex(vtxTrue); | |
326 | Double_t origD[3]; | |
327 | partD->XvYvZv(origD); | |
328 | Short_t charge=partD->Charge(); | |
329 | Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2]; | |
330 | for(Int_t iDau=0; iDau<2; iDau++){ | |
331 | pXdauTrue[iDau]=0.; | |
332 | pYdauTrue[iDau]=0.; | |
333 | pZdauTrue[iDau]=0.; | |
334 | } | |
335 | ||
336 | Int_t nDau=partD->GetNDaughters(); | |
337 | Int_t labelFirstDau = partD->GetDaughter(0); | |
338 | if(nDau==2){ | |
339 | for(Int_t iDau=0; iDau<2; iDau++){ | |
340 | Int_t ind = labelFirstDau+iDau; | |
341 | AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind)); | |
342 | if(!part){ | |
343 | printf("Daughter particle not found in MC array"); | |
344 | return 99999.; | |
345 | } | |
346 | pXdauTrue[iDau]=part->Px(); | |
347 | pYdauTrue[iDau]=part->Py(); | |
348 | pZdauTrue[iDau]=part->Pz(); | |
349 | } | |
350 | }else{ | |
351 | return 99999.; | |
352 | } | |
353 | ||
48efc6ad | 354 | Double_t d0dummy[2]={0.,0.}; |
15917f37 | 355 | AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy); |
356 | return aodDvsMC.ImpParXY(); | |
357 | ||
86d84fa4 | 358 | } |
b176db3b | 359 | //____________________________________________________________________________ |
360 | Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) { | |
361 | // true impact parameter calculation for Dplus | |
362 | ||
363 | if(!partD || !arrayMC || !mcHeader) return 99999.; | |
364 | Int_t code=partD->GetPdgCode(); | |
365 | if(TMath::Abs(code)!=411) return 99999.; | |
366 | ||
367 | Double_t vtxTrue[3]; | |
368 | mcHeader->GetVertex(vtxTrue); | |
369 | Double_t origD[3]; | |
370 | partD->XvYvZv(origD); | |
371 | Short_t charge=partD->Charge(); | |
372 | Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3]; | |
373 | for(Int_t iDau=0; iDau<3; iDau++){ | |
374 | pXdauTrue[iDau]=0.; | |
375 | pYdauTrue[iDau]=0.; | |
376 | pZdauTrue[iDau]=0.; | |
377 | } | |
378 | ||
379 | Int_t nDau=partD->GetNDaughters(); | |
380 | Int_t labelFirstDau = partD->GetDaughter(0); | |
381 | if(nDau==3){ | |
382 | for(Int_t iDau=0; iDau<3; iDau++){ | |
383 | Int_t ind = labelFirstDau+iDau; | |
384 | AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind)); | |
385 | if(!part){ | |
386 | printf("Daughter particle not found in MC array"); | |
387 | return 99999.; | |
388 | } | |
389 | pXdauTrue[iDau]=part->Px(); | |
390 | pYdauTrue[iDau]=part->Py(); | |
391 | pZdauTrue[iDau]=part->Pz(); | |
392 | } | |
393 | }else if(nDau==2){ | |
394 | Int_t theDau=0; | |
395 | for(Int_t iDau=0; iDau<2; iDau++){ | |
396 | Int_t ind = labelFirstDau+iDau; | |
397 | AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind)); | |
398 | if(!part){ | |
399 | printf("Daughter particle not found in MC array"); | |
400 | return 99999.; | |
401 | } | |
402 | Int_t pdgCode=TMath::Abs(part->GetPdgCode()); | |
403 | if(pdgCode==211 || pdgCode==321){ | |
404 | pXdauTrue[theDau]=part->Px(); | |
405 | pYdauTrue[theDau]=part->Py(); | |
406 | pZdauTrue[theDau]=part->Pz(); | |
407 | ++theDau; | |
408 | }else{ | |
409 | Int_t nDauRes=part->GetNDaughters(); | |
410 | if(nDauRes==2){ | |
411 | Int_t labelFirstDauRes = part->GetDaughter(0); | |
412 | for(Int_t iDauRes=0; iDauRes<2; iDauRes++){ | |
413 | Int_t indDR = labelFirstDauRes+iDauRes; | |
414 | AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR)); | |
415 | if(!partDR){ | |
416 | printf("Daughter particle not found in MC array"); | |
417 | return 99999.; | |
418 | } | |
419 | ||
420 | Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode()); | |
421 | if(pdgCodeDR==211 || pdgCodeDR==321){ | |
422 | pXdauTrue[theDau]=partDR->Px(); | |
423 | pYdauTrue[theDau]=partDR->Py(); | |
424 | pZdauTrue[theDau]=partDR->Pz(); | |
425 | ++theDau; | |
426 | } | |
427 | } | |
428 | } | |
429 | } | |
430 | } | |
431 | if(theDau!=3){ | |
432 | printf("Wrong number of decay prongs"); | |
433 | return 99999.; | |
434 | } | |
435 | } | |
436 | ||
437 | Double_t d0dummy[3]={0.,0.,0.}; | |
438 | AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy); | |
439 | return aodDvsMC.ImpParXY(); | |
440 | ||
441 | } | |
442 | ||
443 | ||
444 | ||
686e4710 | 445 | //____________________________________________________________________________ |
446 | Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) { | |
447 | // | |
448 | // Correct the number of accepted tracklets based on a TProfile Hist | |
449 | // | |
450 | // | |
451 | ||
452 | if(TMath::Abs(vtxZ)>10.0){ | |
453 | printf("ERROR: Z vertex out of range for correction of multiplicity\n"); | |
454 | return uncorrectedNacc; | |
455 | } | |
456 | ||
457 | if(!estimatorAvg){ | |
458 | printf("ERROR: Missing TProfile for correction of multiplicity\n"); | |
459 | return uncorrectedNacc; | |
460 | } | |
461 | ||
462 | Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ)); | |
463 | ||
464 | Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1); | |
465 | ||
466 | Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM)); | |
467 | ||
468 | return correctedNacc; | |
469 | } | |
470 |