Correct formula for the error on significance
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliVertexingHFUtils.cxx
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 <TRandom.h>
18 #include <TProfile.h>
19 #include <TClonesArray.h>
20 #include <TH1F.h>
21 #include <TH2F.h>
22 #include <TF1.h>
23 #include "AliAODEvent.h"
24 #include "AliAODMCHeader.h"
25 #include "AliAODMCParticle.h"
26 #include "AliAODRecoDecayHF.h"
27 #include "AliVertexingHFUtils.h"
28
29 /* $Id$ */
30
31 ///////////////////////////////////////////////////////////////////
32 //                                                               //
33 // Class with functions useful for different D2H analyses        //
34 // - event plane resolution                                      //
35 // - <pt> calculation with side band subtraction                 //
36 // - tracklet multiplicity calculation                            //
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
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 }
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 }
196 //______________________________________________________________________
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 //______________________________________________________________________
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
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
354   Double_t d0dummy[2]={0.,0.};
355   AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
356   return aodDvsMC.ImpParXY();
357
358 }
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
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