Correct formula for the error on significance
[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>
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
41ClassImp(AliVertexingHFUtils)
42
43//______________________________________________________________________
44AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
45 fK(1),
46 fSubRes(1.),
47 fMinEtaForTracklets(-1.),
48 fMaxEtaForTracklets(1.)
49{
50 // Default contructor
51}
52
53
54//______________________________________________________________________
55AliVertexingHFUtils::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//______________________________________________________________________
67void 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//______________________________________________________________________
85Double_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//______________________________________________________________________
99Double_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//______________________________________________________________________
105Double_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//______________________________________________________________________
148Double_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//______________________________________________________________________
161Double_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//______________________________________________________________________
168Double_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//______________________________________________________________________
176Double_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//______________________________________________________________________
184Int_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 197Int_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//______________________________________________________________________
210Int_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//______________________________________________________________________
225Int_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 240void 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//____________________________________________________________________________
317Double_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//____________________________________________________________________________
360Double_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//____________________________________________________________________________
446Double_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